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Abstract 


A class of explicit multistage time-stepping schemes tvith centered 
spatial differencing and multigrid is considered for the compressible 
Euler and Navier-Stokes equations . These schemes are the basis for 
a family of computer programs (flow codes with multi grid (FLO MG) 
series) currently used to solve a wide range of fluid dynamics prob- 
lems , including internal a nd external flows. In this paper, the com- 
ponents of these multistage time-stepping schemes are defined , dis- 
cussed, and in many cases analyzed to provide additional insight into 
their behavior . Special emphasis is given to numerical dissipation, 
stability of Runge-Kutta schemes, and the convergence-acceleration 
techniques of multigrid and implicit residual smoothing . Both the 
Baldwin and Lomax algebraic equilibrium model and the Johnson 
and King one-half equation nonequilibrium model are used to estab- 
lish turbulence closure Implementation of these models is described. 


1. Introduction 

Computational fluid dynamics (CFD) is a multidisciplinary field involving fluid mechanics, 
numerical analysis, and computer science. The evolution of CFD over the last three decades has 
fostered a broad range of methods for computing the aerodynamics of flight vehicles. At cruise 
flight conditions, a variety of approximate techniques are applied by the aircraft industry when 
designing flight vehicles. 

With in vise id and irrotational flow assumptions, versatile and reliable pane] methods and 
nonlinear potential equation solvers are used for aircraft design. To determine viscous effects, 
either an integral or finite -difference approach is employed to solve the boundary -layer equa- 
tions. When the interaction between the viscous and in viscid flow regions is important, the 
computational procedures for these regions are coupled in either the direct mode (i.e., surface 
pressure is specified) or the inverse mode (i.e., surface shear stress in the case of a solid wall is 
specified). Although these computational techniques are efficient and usually provide reasonable 
estimates of viscous effects, they can be difficult to implement for three-dimensional (3-D) flows 
when strong viscous-inviscid interactions occur (such as aircraft wing and body juncture flow). 

In the past few years, substantial improvements were made on the mathematical models of 
aerodynamic prediction techniques used for aircraft design. The Euler equations allow rotat ional 
effects (i.e., vortical structures) and nonisent ropic shock waves and thus provide a better inviscid 
model for flows over aerodynamic configurations. The Navier-Stokes equations model weak and 
strong interact ions between viscous and inviscid flow regions without special consideration. Both 
the Euler and the time-averaged Navier-Stokes equations are currently being introduced into t he 
aircraft design process. 

Progress in aircraft design can be attributed to several factors. A primary factor is the 
considerable improvement in the accuracy and efficiency of numerical algorithms used to solve the 
Euler and Navier-Stokes equat ions. Another factor is the significant advancements in computer 
memory capacity and processing times. Although new technologies in computers and computer 
science will continue to help decrease processing times, the need still exists for strong effort to 
increase the robustness, accuracy, and efficiency of the flow solvers to allow their use in analysis 
of complex fluid dynamics phenomena and aircraft design. 

An extensive range of numerical algorithms was developed during the last decade to solve 
the Euler and Navier-Stokes equations. These numerical algorithms can be classified by 
the type of time-stepping scheme and the type of spatial- discretization scheme used. Both 



explicit and implicit time-stepping schemes have been constructed. Explicit schemes require 
less computational storage and a lower number of operations for time integration, but have a 
stricter limit on the allowable time step. If temporal and spatial differencing are decoupled, 
both schemes are amenable to a variety of convergence-acceleration techniques for steady-state 
problems. The explicit multistage Runge-Kutta scheme of Jameson, Schmidt, and Turkel (ref. 1) 
and the implicit approximate factorization (AF) scheme of Beam and Warming (ref. 2) are two 
schemes that employ temporal and spatial decoupling. The multistage schemes, in conjunction 
with local time stepping and other convergence enhancements (ref. 3), and the AF scheme, with 
local time stepping and diagonalization of the implicit operator (ref. 4), are efficient schemes for 
the Euler equations. 

Central and one-sided differencing have been considered for the spatial derivatives in the flow 
equations. When selecting one type of differencing over another, it is important to understand 
the dominating design criterion for central and upwind schemes. When constructing a central 
difference scheme, the principal underlying guideline is to minimize the arithmetic operation 
count while simultaneously maintaining the highest possible accuracy. The multistage schemes 
and Lax-Wendroff schemes (refs. 5-11) are currently the most widely used explicit algorithms 
with central spatial differencing. The AF scheme is the most frequently used implicit scheme 
with centered differencing. 

A primary objective of an upwind scheme is to capture flow discontinuities such as shock 
waves using the minimum number of mesh cells. To accomplish this, many upwind schemes 
utilize the signs of the slopes of characteristics to determine the direction of propagation of 
information, and thus, the type of differencing for approximating spatial derivatives. Two 
procedures for constructing upwind schemes for hyperbolic systems of conservation laws are 
the flux vector splitting scheme of Van Leer (ref. 12) and the flux difference splitting scheme of 
Roe (ref. 13). Upwind schemes have become popular because of their shock -capturing capability. 
Generally, upwind schemes represent shock waves with two interior cells rather than the three or 
four interior cells usually needed by central difference schemes. However, upwind schemes can 
require as much as twice the computational effort. 

Multistage time-stepping schemes with central differencing for spatial discretization on both 
structured and unstructured meshes are now being used to solve the Euler equat ions for flows over 
complex configurations, including airplanes (refs. 14 and 15). Members of this class of algorithm 
have also been extended to allow the solution of the compressible Navier- Stokes equations in 
both two and three dimensions (refs. 16 and 17). Including convergence-acceleration methods, 
such as local time stepping and constant coefficient implicit residual smoothing (which extends 
the explicit time step limit), has made these solvers reasonably effective. Significant performance 
improvements are achieved principally by using the multistage scheme as a driver of a multigrid 
method. The multigrid method involves a sequence of successively coarser meshes and enhances 
the convergence rate and the robustness of the single-grid scheme. In reference 18, a three- 
stage Runge-Kutta scheme with multigrid was successfully applied to the two-dimensional (2-D) 
Navier-Stokes equations. Then, both Swanson and Turkel (ref. 19) and Martinelli and Jameson 
(ref. 20) demonstrated that the type of convergence behavior described in reference 18 could be 
substantially improved. The multigrid procedure was used to solve flow over a wing (ref. 21). 
Significant performance improvements detailed in reference 21 were obtained (refs. 22 and 23) 
by closely following and extending the ideas developed in the 2-D solvers (refs. 19, 20, and 24). 

This paper describes an efficient and versatile class of central difference, finite-volume 
multigrid schemes for the 2-D compressible Euler and Navier-Stokes equations. The elements of 
these schemes are the basis for a family of computer codes (flow codes with multigrid (FLOMG) 
series) developed by the authors that are now being used in both industry and universities. 
These computer codes have been applied to numerous fluid dynamics problems over the last 
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several years, and have been employed as an analysis code in airfoil design procedures (ref. 25). 
The primary purpose of this paper is to discuss, and in many instances, analyze, the components 
of the schemes in these codes. 

Sections 2 and 3 of this report give the flow equations and describe the finite- volume 
formulation for spatial discretization. Three alternatives for numerical approximation of viscous 
stress and heat flux terms are discussed, and the influence of grid stretching on numerical 
accuracy is determined. 

Section 4 of this report discusses artificial dissipation. After outlining the historical devel- 
opment of a form for dissipation, the scalar dissipation model frequently used with the present 
schemes is given in section 4.2. The selection of boundary-point difference operators is an impor- 
tant consideration for the dissipation model. Suitable operators are given, and how local mode 
analysis can often provide an evaluation for a proposed boundary-point difference operator is 
shown. Analysis based on considering the dissipative character of t he discrete system of equa- 
tions is also performed. Section 4.5 examines the intimate connection between the formulation 
of an upwind scheme and a central difference scheme, and a foundation for a matrix dissipation 
model is established. Section 4.6 describes the matrix dissipation model used with the present 
schemes. 

Section 5 discusses the discrete boundary conditions. Section 6 defines t he class of explicit 
multistage time-stepping schemes considered and summarizes their properties. Next, the 
stability of Runge-Kutta schemes for systems of equations is examined. Subsequently, a time- 
step estimate for pseudotime integration of the flow equations is given. Since the temporal 
and spatial discretizations are decoupled, these explicit schemes are amenable to convergence- 
acceleration techniques. Section 7 addresses techniques used in this report, including local time 
stepping, implicit residual smoothing, and multigrid. The initial part of section 7 indicates how 
the discrete system of flow equations is preconditioned with local time stepping. Section 7.2 first 
discusses constant coefficients for implicit residual smoothing, and also presents basic properties 
of residual smoothing. Next., a form for variable coefficients for implicit residual smoothing 
based on stability analysis of a 2-D linear wave equation is introduced. From this form, two 
different formulas for variable smoothing coefficients evolve, and these formulas are compared. 
These variable smoothing coefficients still generally require a time-step estimate that depends 
oil a diffusion limit. The last part of section 7.2 shows how to use variable smoothing coefficients 
to construct new coefficients that can allow removal of the diffusion restriction. Section 7.3 
describes the basic elements of multigrid methods and delineates the salient features of the 
present multigrid algorithm. 

Section 8 discusses turbulence modeling. Both an algebraic equilibrium model and a half- 
equation non equilibrium model are considered. Details for implementation of the turbulence 
models are given. Section 9 states concluding remarks. 

2. Mathematical Formulation 

In this section the integral form of the full Navier-Stokes equations is defined. Boundary 
conditions for the infinite domain problem are then given to complete the general mathematical 
formulation. Section 3 discusses the discrete analogue of the full Navier-Stokes equations, 
section 4 introduces the reduced form of these equations that is frequently solved in aerodynamic 
applications, and section 5 gives boundary conditions for the truncated (finite) domain problem. 

2.1. Equations 

Let p denote the density, u and v represent velocity components in the x and y Cartesian 
directions, respectively; /> is pressure, T is temperature, E is specific total internal energy, and 
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H is specific total enthalpy. If body forces and heat sources are neglected, the 2-D, unsteady 
Navier-Stokes equations can be written in conservative form in a Cartesian coordinate system 
as 


wJl wdv+ J/" dS = ° 


( 2 . 1 . 1 ) 


where t is time, V is the region being considered, 


W = 


T - 


P 

pu 

p V 

IpEj 

puq + jxij. ir-e, 
pvq + pe y + r e y 
pHq + t q - Q J 


and 


q = tie^ + t;e y 


r = (T x e. x e x + T xy e x e y + r yx e y e x + <r y e y e y 
( du 8v\ du 

J ’ = - X \d7 + Ty)- 


T *y - T y* 


du dv \ 
dy dx ) 


du dv 


dv 


(Ty -~ x {d7 + W~ 2 ^ 

Q=t VT= t (ge I+ ^e„) 

E = e + ^(ii 2 + v 2 ) 


H = E + — 

P 

Here, e x and e y are unit vectors of the Cartesian coordinate system (x;y), and n is an outward- 
pointing unit vector normal to the curve S enclosing the region V. Air is the working fluid used 
in this paper. The air is assumed to be thermally and calorically perfect. The equation of state 
is 

p = pRT (2.1.2) 

where R = c p — c iM and the specific heats c y and c v are constant. The quantities // and A 
are the first and second coefficients of viscosity, respectively, and A is taken to be — | p (Stokes 
hypothesis). Either a simple power law or Sutherland’s law can be used to determine the 
molecular viscosity coefficient //. The coefficient of thermal conductivity k is evaluated using 
the constant Prandtl number assumption. The effect of turbulence is accounted for by using the 
eddy-viscosity hypothesis. (See section 8 on turbulence modeling.) 
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2.2. Physical Boundary Conditions 

In the continuum case, either an external or an internal flow problem defined for an infinite 
domain is considered. Thus, appropriate conditions at wall boundaries, which are assumed to 
be solid, must be defined. Later, in the discrete case, finite domains are defined. Suitable inflow 
and outflow boundary conditions mast then be defined. 

For inviscid flows, the tangency (or non penetration) condition 


q n=0 


( 2 . 2 . 1 ) 


mast be satisfied, where q is the velocity vector and n is the unit vector normal to the surface. 
Now, consider the vector momentum equation 


P 


£q 

Dt 


-V;> 


( 2 . 2 . 2 ) 


where Dq/ Dt denotes the substantial derivative of q, and V is the gradient operator. Clearly, 
the substantial derivative of q • n mast vanish along the surface boundary. Therefore 


P 


jUqVj(qn) = 0 


(2.2.3) 


If the inner product of the unit normal and equation (2.2.2) is subtracted from equation (2.2.3), 
then 

pc | (q V) n = n Vp (2.2.4) 

Now, consider the transformation (x, y) — (£ , ?/ ). take i j(x,y) = constant to coincide with 
the surface boundary, and note that the contravariant velocity component \ — — [yc/.J ')(/ + 
(•rp./ -1 ) v (where the subscripts mean differentiation and J is the transformation Jacobian) is 
zero because of equation (2.2.1). Then, from equation (2.2.4) 


P'l 


(*i + y i) 


[(«{*»/ + ;<//,) n + [ y>i u - j y ) ( ( >vj a - puya )] 


(2.2.b) 


For viscous flows, the nonpenetration condition (eq. (2.2.1)) and the no-slip condition 

q t = 0 (2.2.6) 

(where t is t lie unit vector tangent to the surface) must be satisfied. In addition, a boundary 
condit ion is required to determine the surface temperature. For this boundary condition either 
the wall temperature is set to a specified value or the adiabatic condition 

Q - n= 0 (2.2.7) 

is imposed, where Q is the heat flux vector given in equation (2.1.1). 

3. Spatial Discretization 

A finite- volume approach is applied to discretize the equations of motion. The computational 
domain is divided into quadrilateral cells that are fixed in time. For each cell, the governing 
equations can be nondiinensionalized and written in integral form as follows: 

4- // W dx dy + f (F dy-G dx) = [ (F r dy-G v dx) (3.1) 

JJn Jem R< Jem 
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where £2 is a generic cell (or cell area) with <9£2 the cell boundary. In the scaling factor for the 
viscous terms on the right side ofequat.ion (3.1), the quantities 7 , M , and Re are the specific heat 
ratio, Mach number, and Reynolds number, respectively, with M and Re defined by nominal 
conditions. These factors arise because of the choice of nondime nsionalization of the equations. 
The flux vectors are defined by 

r p u 1 


F = 


pu 2 + p 
puv 


puH 


G = 


pv 

puv 

pv' 2 + p 
L pvH 



0 

- 


<*r 


F, = 

Try 



ua d . + VTj.y 



0 

-| 

G( = 

fyx 



1 

+ 

-kZ- 

Oy J 


The independent variables x, y, and / are nondimensionalized as 


x — 


x 

l ref 



/ = i 

tef 


where £t re r = \/p ref/pref? an d the tilde in this section represents a dimensional variable. 
Examples of a reference length are the chord for an airfoil and the throat height for a nozzle 
flow. The thermodynamic variables p, p , and T and the transport coefficients // and k are 
nondimensionalized by their corresponding quantity evaluated at some reference condition. The 
velocity components are scaled by u re j-, and the total quantities E and H are scaled by wjjf f . For 
external flows, the nominal conditions are based on free-stream values, and for internal flows, 
the nominal conditions are based on stagnation values. 

Partition the computational region with quadrilaterals and apply equation (3.1) to each 
quadrilateral. This process is equivalent to performing a mass, momentum, and energy balance 
011 each cell. A system of ordinary differential equations is obtained by decoupling the temporal 
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and spatial terms. In particular, consider an arbitrary quadrilateral (fig. 1, ABCD), and 
approximate the line integrals of equation (3.1) with the midpoint rule. Let the indices (ij) 
identify a cell. Then, by taking W/j as the cell-averaged solution vector, equation (3.1) can be 
written in semidiscrete form as 


w ;j) + £W = 0 (3.2) 

where Q: ,* is the area of the cell, and C is a spatial discretization operator defined by 
C — Cc + £d + £ad* with the subscripts C, D , and AD referring to convection, diffusion, 
and artificial dissipation, respectively. The convective fluxes at the cell faces are obtained by an 
averaging process. The convective flux balance is computed by summing over the cell faces as 

4 

A:' w u = E ( ^V s / (3-3) 

/=! 


with the flux tensor associated with convection given by ( ')/ — F/e^-b G/e^, and for each cell 
face /. the directed length S/ is expressed as 

S / = (A y) { e, d . - (Ax) l e y (3.4) 

where the proper signs of ( A.r)^ and (A j/)y produce an outward normal to the cell face. The 
augmented, convective flux tensor is evaluated as 

(Fch = j(W~q-+ W+q + ), + />, (3.5) 



Figure 1. Finite-volume discretization. 
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where 


— [ 0 (pavg)/©ar (Pavg)/ e t/ ( (/ x l)avg)/ ] 

(Pavg)/ = ^(P~ +P + ) ( 

((pq)avg)/ = j(p“q“ + p + q + )/ 

and the superscripts minus ( — ) and plus (+) indicate quantities taken from the two cell centers 
adjacent to the edge /. The symbol q denotes the velocity vector. In this section, the subscript 
avg always refers to the simple average defined in equation (3.5) for a given edge /. 

The contribution of the diffusive fluxes in equation (3.2) is evaluated as 

4 

C D Wij = ^2(T D ) r S[ (3.6) 

/=! 


where S / is given by equation (3.4), and (^Fd)i ~ (F v )j e a . + (G v )^ey, First-order spatial 
derivatives are in the flux vectors F v and G v . In the present finite-volume method, these 
derivatives are determined using Green's theorem. For example, consider the cell face BC in 
figure 1. The contributions u x and u y to the viscous flux across BC are approximated by their 
mean values as follows: 


( u J-)i x j + 1/2 - ( u x) BC ~ mil, u x dx dy 


= &L'“ ,s 


(3.7) 




ji+ 1/2 


= Mbc= wJJ a u * d * dy 

= -&L" ix 


(3.8) 


where Q f is the area of an appropriate auxiliary cell. 

Three alternatives for computing the diffusion-type terms have been considered. The first 
two approximations for a diffusive flux are obtained with finite-volume methods, and the third 
approximation is determined with a frequently used method based on a local transformation of 
coordinates. A comparison is now’ made between these three choices. 

For the comparison, consider the molecular transport processes associated with cell face BC 
for the x- moment um equation only. Let <& BC = [(Fd)bg ’ s BC']-2 ~ Ay- r yx Ax) BC . In 
one finite- volume method the integration path A 9 B*C 9 D 9 (fig. 1) used in references 16 and 26 is 
considered. Applying the midpoint rule for the required line integrals results in 


®BC 


u <’ J+l + 02 u iJ + 03 U B + 04 u c\ 

. B avg r 1 

+ -qT + (pQVjj + <p-VB + OsrcJ 


(3.9) 
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where 


<P 1 - ~Ay B c Ay pc' + Ax B( • \x B ! C < 

4 

<f >-2 = ^ Ay B c + Ax BC Ax ^ A , 

4 

03 = J Ay BC Ay a i b i + Axb('Ax a i b i 
4 

<P 4 = jAy BC A y ( ’i I y + Ax Bc Ax c , D t 

<pr, - -Ay BC Ax b , c * - Ax BC Ay B i ( ., 

H = | Ay BC Axp^ - A x BC Ayjy A , 

2 

07 = -AyBC Ax A i B > - Axbc Ay A i B i 
og = -Aysc Ax (A iy ~ Ax BC Ay r io' 

1 / x 

= J («/j + «i+l j + «ij + l + «i+l J+l ) 

l / x 

= - (Uj j + Ui-I j + Uij + 1 + M.;-1 j+i ) 

with r/y and v(< defined similarly, and 

Axbc - X C - X B 
Ay BC = Vc- VB 

1 / X 

^avg — 2 v'i'J + Vi j' + l) 

^ = 2 (^' j J+i) 

Martinelli (ref. 27) introduced a different integration path for calculating the viscous terms 
(delineated as DFCE in fig. 2). Integrating around the boundary BFCE with the trapezoidal 
rule results in 


$BC = 


£ avg 

2ft" 


~Ay BC + A-Cflc J Kj+1 - «ij) - ( jAxbc Ay BC j {v;j + 1 - r; 


+ 


+ 


havg 

2ft" 

havg 

2ft" 


-AyyF Ay B c 


("/( - «r) 


—Ax£jr Ay B c ~ A y FF Axbc) (>’B - if) 


(3.10) 
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Figure 2. Alternative integration path for physical diffusive fluxes. 


where 

Axef= Xij+I-Xfj 

4, vef = Vij+i - in,] 

and fi” is the area of the region enclosed by BFCE. The area Q." is given by 

~ 2 (^ X CB 4 Vef - & X EF &Vcb) 

All other quantities in equation (3.10) are defined the same as in equation (3.9). The form of 
®BC §i ven h y equation (3.10) is much more compact, requiring fewer arithmetic operations than 
the form of <I> n( ■ given by equation (3.9). 

A third approach for computing the diffusion- type terms is based on a local transformation 
from Cartesian coordinates (x,y) to arbitrary curvilinear coordinates (£, y). Derivatives with 
respect to x and y are expanded according to the chain rule for partial differentiation. The 
resulting relation for <f?BC is as follows: 


$bc 


A'avg 

ft' 

/'avg 


£&VCBV( + Uff + ( 7j A yCB x t - &XBCy(J v r f 


BC 


+ 


ft' 


^ 4 yBCVy + 4 xg(ix,^\ it£ + C-Ay BC x,, - Axg^Vy) l ’£ 


1BC 


(3.11) 


where 

X Z = 4 x CB 
2/f = 4 y CB 
**/ = 4 x EF 

V>l = &VEF 
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With a uniformly spaced computational domain (A£ — A?/ = 1), &bc m equation (3.11) is the 
same as <J> in equation (3.10), except for the area factor. For a Cartesian mesh, the expressions 
for $BC in equations (3.9), (3.10), and (3.11) are equivalent. If the streamwise-like differences 
associated with the viscous flux quantities are neglected, which is the thin- layer Navier-Stokes 
assumption, only the terms inside the first set of brackets are retained. 

Note that with the thin-layer formulation, there are viscous contributions to the fluxes at 
faces BC and DA only. The following vector approximates the integrand of the right side of 
equation (3.1) at cell faces BC and DA : 

0 
r l 
T 'i 

w avgTL + t'avg 7 ^ + Q 

where 

T \ = (0iU fl - <t> 2 Vr } ) 

11 avg 

r 2 = (H v n - 

i/iiVg 




1 

< l >-2 = 

4 o 2 

4>:\ — + ;v| 

<f>4 = *1 + yl 


Unless otherwise indicated in the text, the thin-layer form of the equations is solved. 

Significant differences in the numerical solutions have not been observed when applying the 
three methods for approximating the diffusive terms. Notable differences in the numerical 
solutions were not expected when solving the Navier-Stokes (thin-layer or full) equations on 
sufficiently smooth meshes (i.e., meshes without kinks or sudden jumps in mesh intervals). In 
this paper the integration path of Martinelli (ref. 27) is used in the finite-volume method for 
computing the viscous fluxes. With this choice of integration path (ref. 27), the mean values of 
the viscous stresses for a given cell edge are obtained at the midpoint of the edge, even when 
there is a kink in the grid. This is not true for the path used in equation (3.9). Also, with the 
integration path of equation (3.10), there are fewer arithmetic operations required than with the 
path of equation (3.9). 

Theoretical est imates of the order of accuracy of the cell-averaged scheme are now introduced 
based on one- dimensional ( 1-D) analysis using day lor- series expansions. Consider the coordinate 
grid around the location denoted by the index i (fig. 3). Let <? be a test function. The numerical 
values of the first and second derivatives of this function are then given by 


A x-\- -j- A.r_ 
A.r 


+ p>jr- 
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respectively, where the derivatives in the expansions are evaluated at i. The approximations of 
equations (3.12) and (3.13) are zeroth-order accurate on arbitrarily stretched meshes. However, 
assuming a constant stretching factor of the grid (i.e., / 3 = Ax^+fAx = constant ), the following 
relations are obtained: 

Ax = Ax 4 

13 

Aj _ = —Ax ( 1 + — ^ 

2 V li) 

= Ax j3 

Az+ = ±Ax(l+/?) 

For viscous flows, grids with constant stretching factor (3 are often used. If these grids are refined 
by doubling the number of points, then 

&j = V&< i + 

where l3 > 1 and the subscripts / and c refer to fine and coarse grids, respectively. To estimate 
the error reduction when refining the stretched mesh, the approximation 

/*« 1 + CjAx (3.15) 

Ls used. Then, if the quantities in equations (3 .14) and (3 .15) are substituted into equations (3.12) 
and (3.13), respectively, the result is 

(t'-’Unuui = 4>r + \<!>r {C a Ax) 2 + Ax 2 + 0(Ax 2 ) (3.16) 
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and 


(^xx ) limn “ 


<f>xx + ~<f>xx {Cj Ax) 2 + -C^xxjr A* 2 + —<t>xxxx A a: 2 + 0{ Aa* ,J 


(3.17 


respectively, for Aar << 1. Thus, second-order accuracy is achieved for the in viscid and viscous 
terms in the flow equations on smoothly stretched meshes. Additional discussion of accuracy is 
found in reference 28. 


4. Artificial Dissipation 

The basic finite-volume scheme described in section 3 contains no dissipative terms in the 
case of in viscid flows. To prevent oscillations near shock waves or stagnation points, artificial 
dissipation terms are added to the governing discrete equations. The introduction of appropriate 
dissipation in the vicinity of shock waves satisfies an entropy condition. In gas dynamics, an 
entropy condition can be the second law of thermodynamics, which states that the physical 
entropy cannot decrease. The entropy condition guarantees the uniqueness of weak solutions 
(i.e., solutions containing shock waves) and thus ensures a physically correct solution. (See 
ref. 29 for further discussion.) 

Another type of dissipation term is added to the discrete flow equations. This term is included 
to provide background dissipation, which is important for converging the numerical scheme that 
will be used to compute flow solutions. These dissipation terms also prevent odd-even point 
decoupling (i.e., creation of sawtooth, or plus-minus waves, with wavelength of two times the 
mesh spacing). For viscous flows, dissipative properties are present because of diffusive terms. 
However, because of the nonlinearity of the equations of motion, the physical dissipation may 
not be sufficient to guarantee stability, especially for the highly stretched meshes generally used 
to resolve the steep gradients in shear layers. Thus, artificial dissipation is also included in 
viscous regions to maintain the stability and robustness of the numerical procedure. 

In this section some historical information regarding the form of the artificial (or numerical) 
dissipation model used with many central difference schemes is discussed. This discussion 
describes how the model evolved, and provides a rudimentary understanding of the model. 
Next, the basic dissipation formulation and various modifications that have been investigated 
are discussed. Boundary-point difference stencils are required for the dissipation model. Several 
stencils are considered and analyzed. Next, the intimate connection between the formulation 
for an upwind scheme and a central difference scheme is examined, establishing a foundation for 
a matrix dissipation model. Section 4.G presents the matrix dissipation model currently used. 
This model relies upon characteristic decomposition of a flux vector. 

4.1. Development of Dissipation Form 

To simplify the historical notes in this section, consider the 1-D system of hyperbolic 
equations (cTW fdt) + (dF/dx) = 0, where W and F are three-component state and flux vectors, 
respectively. Let the 1-D domain be partitioned by intervals defined by A.r = Xj+\ /$ ~~ •**/- 1/2’ 
where the indices refer to interface points for adjacent intervals. Suppose the Lax-Wendroff 
scheme is applied as follows: 


where the superscript n 


F '+l/2 


W" + l = W" - T( F,- +1/2 - F;_!/2 ) 
indicates time level, r ~ Af/Ax, and the interface flux is 

- FfiV/2 = j(F'+ F/+1 ) - TT+ i/j(Wffl - W;) 


(4.1.1) 


(4.1.2) 
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with A representing the flux Jacobian matrix (an element Ajj, - dFj/dWj.). All quantities are 
evaluated at time level n unless noted otherwise. 

In the initial work of computing flows with shock waves by using the Lax-Wendroff scheme, 
the solutions contained oscillations in the vicinity of the shock wave. Then, Von Neumann and 
Richtmyer (ref. 30) introduced an additional dissipation term to remove shock wave oscillations. 
Including this term, equation (4.1.2) is rewritten as 


F i+l/2 = F ™/ 2 - ^+,/2 (W ‘+' - 

or in the continuum 

F = F*- A xd { V^- = F* — D (2) 

OX 

where F* is the physical flux function. The term is often called an artificial (or numerical) 
viscosity and plays the role of a control function. Hirsch (ref. 29) showed that the form of D< 2 ) 
considered by Von Neumann and Richtmyer (ref. 30) can be written for a system as 


D<2) = f ( 2 )Ax 2 $ 


dW dW 

dx dx 


(4.1.3) 


where the coefficients > 0 and can depend on the mesh index i. and each element of D ! 
depends on the corresponding element of W. Now, suppose the flux difference is computed by 
— f /— 1/2) ‘ n equation (4.1 .1) using equation (4.1.3). Then, in the case of the continuum, 
the total dissipation is given by 


dJ 2 ? = £ (2) a * 3 — 

lot dx 




dW 

dx 


<9W 

Ox 


(4.1.4) 


This dissipation term can be characterized as third order. However, Ax 1 appears in equa- 
tion (4.1.1), so effectively, equation (4.1.4) defines a second-order term. 

In 1975, MacCormack and Baldwin (ref. 31) appended a dissipation term for shock capturing 
to the 1969 scheme of MacCormack (a two-step Lax-WendrofF type scheme (ref. 32)). This 
dissipation term was introduced to remove oscillations at shock waves caused by the spatial 
differencing of the MacCormack scheme. This dissipation term is proportional to a second 
difference of the pressure and is given by 

D< 2 > = c( 2 >A a - 3 ^i±I 
4p 


d l p 


dx 2 


d W 

dx 


(4.1.5) 


In smooth regions of a flow field, the product of Ax * and the dissipative flux balance d[ 2 ) Is 

third order, while the product of Ar“* and d[ 2 | is first order in the neighborhood of a shock 
wave. 

As indicated, numerical dissipation is not only important in capturing discontinuities, it is 
also generally required to maintain stability and provide necessary background dissipation for 
convergence. In 1976, Beam and Warming (ref. 2) added to the explicit side of their implicit 
approximate factorization (AF) scheme with what they called a fourth-order dissipation term 
to damp high-frequency error components. With the Lax- Wendroff scheme, this fourth-order 
dissipation term would appear as 


d| 4 ! = -c(4) A ^ 

tol dx4 


(4.1.6) 
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It seems more appropriate to define the order of the dissipation relative to the spatial discretiza- 
tion of the physical terms in the flow equations. So when considering Aj ,_1 d[^J, the dissipation 
term is third order. At any rate, a fourth-difference dissipation is included, and along with the 
second-difference term of equation (4.1.5), provides the basic ingredients for constructing a com- 
plete, adaptive dissipation model. The second-difference term of equation (4.1.5) allows shock 
capturing without oscillations, while the linear fourth-difference term of equation (4.1.6) pro- 
vides the important background dissipation. The critical element missing is a switching function 
that would turn on the appropriate dissipation form in a region and turn off the dissipation form 
that is not the desirable type (i.e., near shocks of eq. (4.1.5) should dominate, with D^) 
of eq. (4.1.6) negligible, while in smooth regions. of eq. (4.1.6) should dominate, with 
of eq. (4.1.5) negligible). In section 4.2, the dissipation model that adds a switching function to 
the two basic types of dissipation terms just, discussed is described. 

4.2. Dissipation Model 

To permit a complete description of the dissipation model, the two-dimensional Euler 
equations are now considered. The dissipation is based on the model introduced by Jameson, 
Schmidt, and Turkel (ref. 1) that defined a suitable switching function (at least for transonic 
and low supersonic flow) to allow blending of the second and fourth differences. According to 
the nonlinear model (ref. 1), the quantity j in equation (5.2) is expressed as 

C AD Wij = -(£>? - i)} + D 2 h - dJ)W,j (4.2.1) 

where (£,r/) are arbitrary, curvilinear coordinates, and 


(2) 


DfW,j = V* (A ;+1/2jf } + ' /2 .)A f W ; j 


(4.2.2) 


Dt W; ; = Vc 

£ * J <N 


O;. 


'+ 1/2 J 4 + 1 / 2 , j ^ ^ ^ ^ 


W; 


' J 


(4.2.3) 


where i and j are indices (for a cell center) associated with the £ and ?/ directions, respectively, 
and A^ and are forward and backward difference operators in the £ direction, respectively. 
The definitions are similar in the r; direction. The variable scaling factor A is defined as 


A '+1/2J - 9 W,-; + ( A {), + 1J + + (A '/t+lJ 


(4.2.4) 


where Xc and X f( are the largest eigenvalues in absolute value (i.e., spectral radii) of the flux 
Jacobian matrices associated with the Euler equations. These spectral radii are given by 


-V — \ u V‘i - + C \J ~y #/ + 

> 

\,, = \vj'£ - vyc I + CyJ. J‘i + yj? j 


(4.2 J) ) 
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where u and v are Cartesian velocity components, and c is the speed of sound. The coefficients 
and e use the pressure as a sensor for shocks and stagnation points, respectively, and are 
defined as 


d+l/2 J = kV2) ma x( Vi - 1 j , Vj j ^(4-2 j) 


"ij = 


Pi-U ~ 2 PiJ + Pi + lj 
Pi— IJ T j T Pi+lJ 


(4.2.6) 


,(4) 


= max 


0,(k 


(4) _ ^(2) 


where typical values for the constants /c^ and are in the ranges 1/4 to 1/2 and 1/64 to 1/32, 
respectively. This paper shall refer to equations (4.2.1) and (4.2.6) as the Jameson, Schmidt, 
Turkel ( JST) scheme (or dissipation model), and shall designate v as the JST switch. It should 

be mentioned that in reference 1, the coefficient j — max(^j,z/;+i j). The switching 

function v can be interpreted as a limiter, in the sense that it activates the second-difference 
contribution at extrema and switches off the fourth- difference term. Moreover, at shock waves, 
the dissipation is first order, and a first-order upwind scheme is produced for a scalar equation. 
In smooth regions of the flow field the dissipation is third order. 

Th ns, two different dissipation mechanisms are at work, and the switch determines which 
one is active in any given region. For smooth flows, v is small, and the dissipation terms consist 
of a linear fourth difference that damps the high frequencies the central difference scheme does 
not damp. This dissipation is useful mainly for achieving a steady state and is less important 
for time- dependent, problems. In the neighborhood of large gradients in pressure, u becomes 
large and switches on the second-difference viscosity while simultaneously reducing the fourth- 
difference dissipation. This viscosity is needed mainly to introduce an entropy condition so that 
the correct shock relationships are satisfied and to prevent oscillations near discontinuities. For 
subsonic steady-state flow, this viscosity can be turned off by choosing = 0 . 

The isotropic scaling factor of equation (4.2.4) is generally satisfactory for inviscid flow 
problems when typical inviscid flow meshes (i.e., cell aspect ratio 0(1)) are used. The isotropic 
scaling factor can create too much numerical dissipation in cases of meshes with high-aspect- 
ratio cells. The adverse effect of high-aspect-ratio cells is an important consideration for high 
Reynolds number viscous flows, where a mesh providing appropriate spatial resolution can have 
cell aspect ratios O(IO^). In an effort to improve this cell aspect ratio situation and obtain 
sharper shock resolution on a given grid, Swanson and Turkel (ref. 19) replaced the isotropic 
scaling factor of equation (4.2.4) with the anisotropic scaling factor 

x - 1 

A similar scaling is used in the 77 direction. 

The anisotropic scaling idea was motivated by the scaling of dissipation occurring in 
dimensionally split, upwind schemes (i.e., the flux vector split scheme (ref. 12 ) and the 
approximate Ri email n solver (ref. 13)). Anisotropic scaling is often referred to as individual 
eigenvalue scaling. While the accuracy is improved with equation (4.2.7), particularly with 
respect to shock resolution, individual eigenvalue scaling in the streamwise (£) direction can 
be too severe for a standard multigrid algorithm. Moreover, the effectiveness of the multigrid 
driving scheme in damping high frequencies in the £ direction can be significantly diminished, 
resulting in a much slower convergence rate. 


(A f ) +(A f ) 


(4.2.7; 
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An alternative to the individual eigenvalue scaling was proposed by Martinelli (ref. 27), and 
considered by Swanson and Turkel (ref. 19). This modified scaling factor, which is a function of 
mesh -cell -aspect ratio, is defined as 




(4.2.8) 


where 


=^f lJ f( r ) (\)/j 
<f } ij( r ) = 1 + r /J 


(4.2.9) 


Here, r is the ratio A r/ /A<f, and the exponent £ is generally taken to be between 1/2 and 2/3. 

In the normal direction (?/), (A , ; );j — 1 )( A /; ); j is defined. Thus, the scaling factor 

of equation (4.2.8) is bounded from below by equation (4.2.7), and bounded from above by 
equation (4.2.4). As demonstrated in references 19 and 20, the scaling factor of equation (4.2.8) 
produces a significant improvement in accuracy for high-aspect-ratio meshes, and permits good 
convergence rates with a multigrid method. The scheme in this paper uses this modified scaling 
factor. 


The impact of the dissipation form on the energy of a. system of flow equations is now 
examined. For simplicity, consider the 1-D, time- dependent Euler equations, with numerical 
dissipation terms given by 


where 


r — 1 


<9W <^F 




D f w = £ 

£ X V ) 




f(4) 


yw 

oe 


(4.2.10) 


and A ; r ( ^ = Af^K and J — * is the inverse transformation Jacobian. Form the 

inner products of W r , with T denoting transpose, and both sides of equation (4.2.10), and 
then integrate over a domain Q. After integration by parts and neglecting boundary terms, the 
equation 

7^-11 W|| 2 = flux term + J — 7^j (4.2.11) 


Ls obtained, where 

|| W|| 2 = f W 2 







The second-difference dissipation term /W only decreases the L 2 norm of the solution vector 
(i.e.. it decreases the energy of the system), and thus, is strictly dissipative. The fourth -difference 
dissipation term contains a dispersive part and a dissipative contribution. 
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An alternative form for the third-order dissipation term (the last term in eq. (4.2.10)) is 


Di w 


_£ 2 

di 


, 10 ^' 1 

OF) 


which can be written in the discrete case as 


Ac 


(A,~4V s cAc 


W i 


(4.2.12) 


Tliis modified form produces only dissipative contributions. If Q, = A/ej^V^ A^W/, then 
Vf A^Q, = <*>Q/_|_i/ 2 — ^Q/-l /25 where h is the standard, centered-difference operator. Therefore, 
the form of equation (4.2.12) is conservative, provided A and are evaluated at the cell centers 
rather than at the cell faces. In reference 19, numerical tests were performed with the dissipation 
terms of equations (4.2.3) and (4.2.12). For steady state, there seems to be no consistent benefit 
for either convergence or accuracy when using the form of equation (4.2.12). Based on these 
results, the form of equation (4.2.3) is still used. However, for unsteady flows, equation (4.2.12) 
may offer an advantage because of the absence of dispersive effects that can cause phase errors. 

Until now, a scalar viscosity in which the viscosity is based on differences of the same quantity 
advanced in time has been considered. (See eq. (4.2.10).) The disadvantage is that the total 
enthalpy is no longer constant in the steady state, even when total enthalpy should be identically 
constant for the inviscid equations. The total enthalpy is constant for the steady-state Euler 
equations because the energy equation is a constant multiple of the continuity equation when H 
is constant. Hence, reference 1 suggests that the dissipation for the energy equation be based on 
differences of the total enthalpy rather than the total energy. Thus, a typical situation in one 
dimension is to replace equation (4.2.10) for the energy equation by 


dpE 

~ 


+ 




( 4 . 2 . 13 ) 


where pH = pE-Y p . Reference 33 shows that equation (4.2.13) indeed yields a constant total 
enthalpy, but that the entropy tends to be less accurate than if the dissipation term for the energy 
equation is based on differences of pE rather than pH . Thus, both choices have advantages and 
disadvantages. The total enthalpy formulation is used in this paper. 


4.3. Boundary Treatment of Dissipative Terms 

In a cell-centered, finite-volume method, the first and last cells in each coordinate direction are 
auxiliary cells where the flow equat ions are usually not solved. The solution in these cells is found 
by a combination of the given physical boundary conditions and numerical bound ary conditions. 
Thus, there is no difficulty evaluating the second-difference dissipation term at the first or last 
interior cell in a given coordinate direction. In the case of the fourth- difference dissipation term, 
the treatment must be modified at the boundaries of the physical domain because only one 
layer of auxiliary cells is considered. Moreover, the standard five-point difference stencil must 
be replaced at the first two interior mesh cells relative to a wall boundary; thus, one-sided or 
one-sided biased stencils are used at these cells. The dissipat ive character of these stencils is 
important because it influences both stability and accuracy. For example, if the dissipation is 
too large at a solid boundary, an artificial boundary layer is created in an inviscid flow, and the 
effective Reynolds number for a viscous flow is altered. 


4*3.1. Boundary-point operators . In this section, the two types of discrete boundary-point 
operators (difference stencils) used with the present scheme for solid surfaces are defined. 
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1 j = 1/2 

Figure 4. Boundary- point dissipat ion. 


Next, these operators are evaluated by applying a. local mode analysis. In addition, this 
section shows how this local mode analysis can provide an evaluation of candidate boundary- 
point operators once a basis for comparison is established. A more complete analysis for 
the boundary-point operators is based on the dissipation matrix for the system of difference 
equations approximating the governing flow equations. Sometimes the dissipation matrix can 
be characterized analytically. In general, the eigenvalues of the dissipation matrix must be 
determined. The approach for analyzing the dissipation stencils is discussed. 

Consider the total dissipation resulting from a numerical flux balance for a mesh cell in a 
particular coordinate direction. Let Wj and dj denote a component of the solution vector W 
and the corresponding total dissipation, respectively. The index j indicates the mesh cell being 
considered. Let dj+ 1/2 an d 1/2 represent the dissipative fluxes at the cell interfaces j + 1/2 
and j — 1 / 2, respectively (fig. 4). At a cell interface (for example, j-h 1/2), let ( denote 
the difference between the solution for the adjacent cells (wj + 1 — ivj). For simplicity, assume 
As( 4 ) = 1. Then, for any cell j 

dj - dj+ 1/2 - dj- 1/2 (4.3.1) 

where the dissipative fluxes are 

dj+i /2 = (AuOj+3/2 - 2( Au-')j+i /2 + (Atr)y_i /2 
dj- 1/2 - (Aw) j+ 1/2 - 2(Au>)j_i/2 + (A«;)j_3/2 

Thus 

dj = ( Aw) j+ 3/2 - 3( Aw) j+ 1/2 + 3(A u - ( Au- )j _ 3/ 2 (4-3.2) 

or 

dj = wj+ 2 - 4wj + i -f (ht)j - 4wj_i + « j _ 2 

Consider the first two interior cells adjacent to a solid hound ary (fig. 4 ) . The total dissipation for 
these- cells is denoted by d> and c/ 3 . At j - 2, a value for (Aw)^ »» Lst ,,e determined because 
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(Au>)]/ 2 is undefined. Also, in this formulation of the boundary-point dissipation stencil, no 
functional dependence on uq is desired because w\ is outside the domain. Hence, a value for 
(Aw^m must also be provided. If 



(Att?) 1/2 = (Aw;) 5/2 ) 
(Aw) 3/2 = (Aw ) 5/2 J 

( 4 . 3 . 3 ) 

then equation ( 4 . 3 . 2 ) gives 

d'2 ~ W4 — 2w 3 -f w ( i 

( 4 . 3 . 4 ) 


d% ~ — Aw4 + 5 w 3 — 2 tv 2 

( 4 . 3 . 5 ) 

These boundary stencils are fairly standard and are used for inviscid flow 
alternative form, which reduces the sensitivity to solid-surface, normal mesh 

’ calculations. An 
spacing for viscous 


flow calculations without compromising stability or convergence, is obtained by replacing 
(Aw) w »th (Aw) = 2(Aw)3/2 - (A w)§/2 and leaving (Aw^^ unchanged. This form 
is given by 

d-2 — W4 — 3w3 + 3u'2 - w\ ( 4 . 3 . 6 ) 

d% = wfy — Aw 4 + 6x03 — Aw 2 + w\ ( 4 . 3 . 7 ) 

For turbulent flows, this boundary dissipation formulation (eq. ( 4 . 3 . 7 )) is advantageous when 
the mesh is fine enough to adequately represent the laminar sublayer region of the boundary 
layer (i.e., at least two points are inside the sublayer). For coarse meshes, this treatment of the 
dissipation can be less accurate than the zeroth-order extrapolation of equations (4.3.3). 


4- 2. Local mode analysis . A local mode analysis is now considered to evaluate the 
relative damping behavior of boundary-cell difference operators. For comparison purposes, the 
interior fourth difference is first characterized. Taking a Fourier transform of equation ( 4 . 3 . 2 ) 
yields Zj(0) = 4 (cos 0 — l) 2 , where Zj(0) is the Fourier symbol of the transformed dj , and 6 

is the product of the wave number and the mesh spacing. Then, zj{6) ~ 6 4 for small 6 , and 
zj( 7r) = 16 . The dissipation of long wavelengths is dictated by the behavior of zj(0) at small 
6 , and the dissipation of short wavelengths is governed by zj( tt). As mentioned initially in this 

section, this simple analysis assumes that = 1 . In practice, the coefficient used in 

the evaluation of f or the fourth-difference dissipation affects the behavior of the boundary 
dissipation stencil. The coefficient /^ 4 ) j s chosen such that the highest frequency is highly 
damped according to a stability analysis using the interior-point stencil. This is important for a 
multigrid method and will be discussed in section 7 . 3 . Near a boundary, the dissipation should 
behave in a similar manner. In this dissipation model, the same value of used for interior 
points of the domain is also used near a boundary. 

A general form of the difference stencils at j = 2 , 3 can be written as 

d j = aw j+ 2 - dw j+i + (0 + 7 “ a ) w j ~ l w j~\ 

The associated Fourier symbol is given by 

- j(@) — \0 + 7 — 2q ( 1 + cos 6)} ( 1 — cos 6 ) 

+ * (7 — 0 + 2 a cos 6 ) sin 0 
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For small 6 this Fourier symbol is replaced by 


0 4 

z j(&) — {0 + 7 — 4 a ) — + a — 



+ i(2a — /? + 7 

(4.3.8) 

and at 6 — tt reduces to 

Zj(n) = 2(/i + 7) 

(4.3.9) 

In the case of equation (4.3.4) 

e 4 •> 

Z2(0)*j-i0 3 

(4.3.10) 

for small 0, with z 2 (^) — 4, an 

d for equation (4.3.5) 



-3 ) ~ i()i 

(4.3.11) 

for small 0, with 23(71*) = 12. 

Note that z 2 {0) and z •>,(()) are not. real. 

Thus, there are both 


dissipation and dispersion near the boundary. For the stencil of equation (4.3.6) 

z 2 (4.3.12) 

for small 0 , with ^(tt) — 8, and for the stencil of equation (4.11.7) 

s 3 (0)«0 4 (4.3.13) 

for small with 23(71-) = 16. Comparing equations (4.3.10) and (4.3.12). which correspond 
to the stencils of equations (4.3.4) and (4.3.6), respectively, shows that, both stencils behave 
the same for the long wavelengths, while equation (4.3.12) is twice as dissipative for the short 
wavelengths. At j — 3, the stencil corresponding to equation (4.3.13) is fourth order on the 
long wavelengths, whereas the stencil associated with equation (4.3.11) is only second order. In 
addition, the symbol of equation (4.3.13) is more dissipat ive on the short wavelengths. Thus, 
the improved accuracy and high-frequency damping observed for t he stencils of equations (4.3.6) 
and (4.3.7) in practice is substantiated with this simple analysis. 

The method of combining the simple local mode analysis with the evaluations just considered 
to quickly evaluate candidate dissipation stencils can now be shown. Consider a different set 
of boundary-point stencils. If Aw is taken to represent either the component pu or pv of the 
solution vector W, and the antisymmetry constraint (Aw)^ = (Aii')^ is imposed for viscous 
flows, then equation (4.3.2) gives 


c 1*2 = W 4 — 5 u j 3 + 7u>2 — 3 w i (4.3.14) 

g? 3 = ir 5 — 4 W4 T 6 u'3 — 4w 2 + «’l (4.3.15) 

The Fourier symbols of equation (4.3.14), using equations (4.3.8) and (4.3.9), are 

z 2 (0) « 20 2 (4.3.16) 

for small 0, with z 2 (tt) =16, and the symbols for equation (4.3.15) are the same as given in 
equation (4.3.13). Comparing equation (4.3.16) with equation (4.3.12) shows that the highest 
frequency is damped better with the proposed stencil, but that the proposed stencil is only second 
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order on the long wavelengths, while the stencil of equation ( 4 . 3 . 6 ) is fourth order, indicating 
that better accuracy is obtained with equations (4.3.6) and (4,3.7). The improved accuracy has 
been verified with numerical experiments (i.e., skin-friction solutions for turbulent airfoil flows 
have been computed on 160 by 32 meshes and compared with high- density -mesh results). 

4.4. Matrix Analysis 

The associated dissipation matrix is examined to determine the numerical dissipativity of a 
discrete system of equations, such as equation (3.2). For simplicity, consider the 1-D system 

where w is a discrete solution vector, and D W is a dissipation matrix corresponding to 
fourth-difference terms. Taking the inner product (the transpose ofw) with each side of 
equation (4.4.1), obtain 1/2 dW 2 /dt = w T . If the quadratic form w T D^w is nonpositive 

definite, then the matrix /l ( tl is strictly dissipative. Moreover, the energy of the system is 
nonincreasing. Assume there are boundaries at j = 3/2 and j = jl + 1/2, and assume j = 2 and 
j = jl are the indices for the first, and last interior points, respectively. Apply the boundary point 
stencils of equations (4.3.4) and (4.3.5) at the first two interior cell centers at both boundaries, 
and the standard stencil everywhere else. The resulting dissipation matrix is given by 




■-1 2 

-1 



“ 




2 -5 

4 -1 







-1 4 

-6 4 

-i 






0 -1 

4 -6 

4 -1 





d (4) = 






(4.4.2) 




-1 

4 -6 
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-1 0 






-1 4 

-6 

4 -1 






-1 
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-1 
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and the 

corresponding solut: 

ion vector 

is given by 





w=[. 

u ’2 

w 4 w 5 

• Wjl -2 

Wjl 

-1 w jl f 

Then 










jl - 1 






( w j+i 

— 2 w, -f it 

'j - 1 

) 2 < o 

(4.4.3) 


Thus, D* 4 ) is strictly dissipative. This same result is obtained by Eriksson and Rizzi (ref. 34). 
For a 10 by 10 matrix with the form of equation (4.4.2), Pulliam (ref. 35) obtains two zero 
eigenvalues. Ideally. should have no zero eigenvalues, since zero eigenvalues can possibly 

produce undamped modes that cause instabilities (ref. 35 ). 

Pulliam (ref. 35) recommends applying a stencil with the weights of equation (4.3.5) at the 
first, interior cell, and a standard stencil with the weights of equation ( 4 . 3 . 7 ) at the second interior 
cell. Then 
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and 


■-5 4 -1 

4 -6 4 -1 

-1 4 -6 4 -1 

0-14-64-1 


D (4) = 


- 14 - 64-10 

-1 4 -6 4 -1 

-1 4 -6 4 

-1 4 -5. 


^ ( t/y+i - ‘2wj + v-j-i )" 
y'=3 

- («g - '2w‘> ) 2 - («-/_! - 2wji) i < 0 


(4.4.4) 


(4.4.5) 


Again, the dissipation matrix is strictly dissipative. Moreover, a 10 by 10 matrix with the 
structure of equation (4.4.4) has zero eigenvalues (ref. 35). However, indications are that 
for a cell-centered, finite- volume formulation, this boundary-point treatment of the dissipation 
wit h the weights of equations (4.3.5) and (4.3.7), although appropriate at inflow and out flow 
boundaries, is generally too dissipative at solid boundaries. Thus the stencils of equations (4.3.4) 
and (4.3.5) are preferred at a wall boundary. 

Now consider the stencils with the weights of equations (4.3.6) and (4.3.7). The dissipation 
matrix is given by 
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(4.4.7) 
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From the quadratic form of equation (4.4.7), it does not directly follow that is nonpositive 
definite, which is generally the case with the quadratic form. If the eigenvalues of a 10 by 
10 matrix with the structure of equation (4.4.6) are determined, one is zero and the others 
are negative. Therefore, the matrix D ^ is nonpositive definite. Although there is one zero 
eigenvalue, the present scheme performs well using the boundary-point operators associated 
with equations (4.3.6) and (4.3.7) at solid boundaries when solving viscous flow problems. 


4.5. The Upwind Connection 

Upwind schemes for solving hyperbolic systems of conservation laws (i.e., Euler equations of 
gas dynamics) generally rely upon characteristic theory to determine the direction of propagation 
of information and. thus, the direction required for one-sided differencing approximations of the 
spatial derivatives. With upwind schemes, shock waves can be captured without oscillations. 
Thus, a successful artificial dissipation model for a central difference scheme should imitate an 
upwind scheme in the neighborhood of shocks. The connection between upwind and central 
difference schemes is now T reviewed. 

Consider the 1-D scalar wave equation 


du du 

W + a te = 


o 


with a constant. The first-order upwind scheme can be written as 


u '' + 1 = Uj - u-h 
J J Ax 


Uj + 1 - Uj (a < 0) 
uj — uj _ | (a > 0) 


(4.5.1) 


where all discrete quantities are evaluated at time level n At unless otherwise denoted. The 
scheme of equation (4.5.1) can be rewritten as 


v/' +1 

J 


U J - 


A/ At 

1 - u j- 1) + ^jA^ Uj+l “ ' lUj + : 


(4.5.2) 


Equation (4.5.2) now contains a central difference term and a second-difference dissipation term. 
Now* consider the system 

0u du 

— + A— = 0 (4.5.3) 


dt 


dx 


where u Ls an JV-component vector. The system case can be converted to a scalar system by 
diagonalizing the N by N matrix A with a similarity transformation A — T _1 AT, where the 
columns of T are the right eigenvectors of A. After diagonalizing equation (4.5.3), and applying 
the scheme of equation (4.5.2), the first-order upwind scheme is given by 


/j+1 ^ 

U C =u ^- a 2A7 


At 


u 


j+1 


U j-l) + l- 4 lj^j( u j+l - 2u j + u j-l ) 


(4.5.4) 


w r here 


\A\ = T\A\T~ l 
A = Diag [| Ai| |A.v|] 


(4.5.5) 
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Note that since A has only three distinct eigenvalues, by using the Cayley- II amilton theorem. \A\ 
can he expressed as a quadratic polynomial in A. The generalization to a system of conservation 
laws is as follows: 



dx 


= 0 


with f being an JV- component flux vector, and 


u " +1 = u, 
J J 


At -(i j+ i-f h i)+T At 


‘2 Ax 


2 Ax 


A j+ 1/2 


( u j + l - u j) “ 


A 


j- 1/2 


( u j 


u 


■J - « 


(4.5.6) 


where the Jacobian matrix A = df/du. and |,4| is defined the same as for equation (4.5.4). The 


matrix 


p4j + l/2| can be computed as either an arithmetic average or a Roe average (ref. 13). 
For transonic, steady flows the differences are negligible and the simpler arithmetic average is 
used. Yee (ref. 36) found that the Roe average yields better results for hypersonic flows. The 
Roe average also seems to give slightly better results for time-dependent problems. 


4.6. Matrix Dissipation Model 

As indicated in sect ion 4.5, high resolut ion of shock waves without oscillations can be achieved 
by closely imitating an upwind scheme in the neighborhood of a shock wave. A key feature of 
upwind schemes is a matrix evaluation of the numerical dissipation. W ith this mat rix evaluation, 
the dissipative terms of each discrete equation (associated with a given coordinate direction) are 
scaled by the appropriate eigenvalues of the flux Jacobian matrix rather than by the spectral 
radius, as in the .1ST scheme. Such a matrix dissipation also allows high resolution of wall 
bounded shear layers (ref. 37). The modifications of the .1ST dissipation model required to 
produce the matrix dissipation model currently used are now presented. 

Consider the two-dimensional, time-dependent Euler equations in the form 


fl(7 _1 W) dF dG 
dt dij 


(4.6.1) 


where F and G are flux vectors, W is the solution vector, and (£. ?/) are arbitrary curvilinear 
coordinates. Define A and B as the flux Jacobian matrices AF/AW and AG /AW. respectively. 
By extending the scheme given in equation (4.5.6) to two dimensions, it follows that t he mat rices 
|A| and | B\ must be the scaling factors in a matrix dissipation model. Now, consider the .1ST 
dissipation model. The necessary modification to the contributions for the £, direction of the 
artificial dissipation term defined by equation (4.2.1) is to substitute matrix |A| for the eigenvalue 
scaling factor A in equations (4.2.2) and (4.2.3). For the ij direction, £ and matrix |.4 1 are 
replaced by i] and matrix |fl|, respectively. Next, define explicitly the form for the matrix |/1|. 
Let A = Diag [Ai A2 A3 A3] with 


A 1=9 + \J « i + “2 c 

A 2 = 9 - \f a { + a -i r 
A3 = <i 
« 1 = J ~ 1 6- 

<j — a [ ti + a-> v 
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Then, 


HI = |A 3 |/ + (Hi!±JM_| A3 |) (V^ £l+ 1 £.; 
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(4.6.2) 


Here, H is the total enthalpy, and <j> — (lA + r;^)/2. Note that for the matrices E j , each row 
is a scalar multiple of the other rows (except for zero rows). Hence, to find the product EfW. 
simply find one element of the product Ej W, and the other rows are then scalar multiples of 
that element. Because of the special form of matrix \A\ for any A], A 2 , and A 3 , an arbitrary 


vector x can be multiplied by matrix |/1| very quickly. That Is, calculate 


V+l/2 


(Wj+i-Wj 


directly rather than calculate 


Ah- 1/2 


computed the same way as matrix \A\ 


j and multiply a matrix by a vector, 
by simply replacing £ with i ]. 


The matrix \B\ is 


26 


In practice, Ai,A* 2 , and A 3 cannot be chosen as given above. Near stagnation points, A 3 
approaches zero, while Aj or A 2 approach zero near sonic lines. A zero artificial viscosity creates 
numerical difficulties. Hence, these values are limited as 

I Ail = max[|A! |, V n p{A)} ' 

P(A) = \q\ + cyla{ + al 

> 

IA 2 I = max[|A 2 |, Vnp(.4)j 

IA 3 I = max[|A 3 |, l}p(A )] j 

where the linear eigenvalue A 3 can be limited differently than the nonlinear eigenvalues. The 
parameters Vn and V( were determined numerically. Various values were evaluated by comparing 
their corresponding computed solutions based on the sharpness of shock weaves captured (without 
producing oscillations) and convergence rate of numerical scheme. Based on this evaluation, a 
good choice for V n and 1/ is 0.2. However, in reference 37, accurate coarse-grid solutions for a 
low-speed, high Reynolds number (5x 10'^) laminar flow over a flat plate were not obtained with 
V/ = 0 . 2 . Accurate coarse-grid results (i.e., 5 to 10 points in boundary layer) w 7 ere computed with 
V( — 0.01 for the direction normal to the plate, and l ) — 0.2 for the streamwise-like direction. 

Thus far, A /+j /2 j in equations (4.2.2) and (4.2.3) has been replaced by a matrix while leaving 

the limiters and as scalars. Also, and can be introduced into the diagonal matrix 
A, allowing different limiters to be chosen for different characteristic variables. For example, 
the limiter may be based on pressure for the nonlinear waves. However, the pressure is smooth 
through a contact discontinuity. Hence, aswutch based on temperature maybe more appropriate 
for the linear wave. Different ntesh scalings, and thus different <j>{r) for the linear and nonlinear 
waves, could also be used. 

5. Discrete Boundary Conditions 

An important element when developing an accurate and efficient algorithm for solving the 
Euler and Navier-Stokes equations is selection of proper boundary conditions. The choice of 
conditions must be consistent w ith physical const raints of the problem of interest and the interior 
discrete formulation. Moreover, the physical conditions generally must be supplemented with a 
sufficient number of numerical relations to allow determination of all dependent variables. 

In addition to defining the conditions at solid or porous w 7 all boundaries, the infinite domain 
problem must be adequately simulated for external airflows. External airflow 7 simulation is usu- 
ally done by delineating boundaries at some distance from the primary region of consideration, 
and then prescribing suitable boundary conditions for that location. In the case of a lifting 
airfoil, the outer boundary position must be far enough away from the airfoil not to compromise 
the development of the lift. For example, 5 airfoil chords w T ould be too close, whereas 20 chords 
would be satisfactory if the far-field vortex effect (ref. 38) is considered. Even for in viscid, non- 
lifting airflow over a circular cylinder, an outer boundary placed too close to the cylinder can 
cause inaccurate prediction of the airflow 7 over the aft portion of the cylinder. 

At a solid boundary, a row of auxiliary cells is created exterior to the domain of the airflow . 
By approximating the normal pressure gradient of equation (2.2.4) with a three-point centered 
difference at the surface, the auxiliary cell pressure is obtained. The density at this cell is 
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equated to the density at the first point off the surface. The tangency condition is enforced by 
determining the Cartesian velocity components from 


u 


> 

I 


Hi 

V 

U 



IV 

m Qn_ 


where £ is the coordinate aligned with the surface boundary, qt and q n are the tangential and 
normal velocity components, respectively, the subscript w means wall, and the indices (i, 1) and 
(b 2) refer to the centers of the auxiliary and the first interior cells, respectively. The overbar 

means the quantity is divided by y/(ar| + y|). Finally, the total internal energy is computed 
using the relation 

1 loo 

pE = -p + -p(u z + v l ) 

7 — 1 / 

In the case of viscous flows, the no-slip condition is required, and is imposed by treating the 
Cartesian velocity components as antisymmetric functions with respect to the solid surface. 
Thus 

u i,l = -«/, 2 

Vi, 1 = -Vj/2 


The surface values of pressure (p) and temperature (T) are computed using the reduced normal 
momentum and energy equations 



dT 

drj 


= 0 


> 


( 5 . 1 ) 


where 7 / is the coordinate normal to the surface. As part of the boundary conditions, the option 
to specify the wall temperature instead of imposing the adiabatic condition of equations (5.1) is 
included. 


To compute the unknown flow variables at the outer boundary of an external aerodynamics 
problem, characteristic theory, some simplify ing assumptions, and the concept of a point vortex 
are used. In appendix A, a point on the outer boundary and the two-dimensional Euler equations 
are considered. Then, assuming a locally homentropic flow, the one-dimensional equations of 
gas dynamics are derived (for completeness) for the direction normal to the boundary. The 
elements of the solution vector are proportional to the local tangential velocity component and 
the Riemann invariants 


R + — q n + 


2c 

7 - 1 


and 


R — <7,i 


2c 

7 - 1 


respectively, where the tangential and normal velocity components are defined as 


\A*f + 


+ y{) 
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respectively. This set of dependent variables, the honientropic assumption, and characteristic 
theory are used to determine the unknown flow variables. 

To compute the discrete solution at the outer boundary points (as for the wall boundary), a 
row of auxiliary (boundary) cells exterior to the domain is introduced. Then, at a boundary cell, 
the normal velocity component q n and the speed of sound c are computed from the relations 

% = i(/? + + /n 

c = l^-(n+-ir) 

where the characteristic variables /?+ and /?“ are appropriately determined. Assume that 
the flow normal to the boundary is subcritical. If inflow occurs, the characteristic variables 
corresponding to the ingoing characteristics are specified. Since this is actually a two-dimensional 
system, an additional quantity must be given. It follows directly that the entropy s should be 
specified (the flow is assumed to be locally homentropic). In practice, for convenience define 
s* — pj p~ ! , which has the same functional dependence as entropy, and use this variable in place 
of entropy. So, for an inflow situation, set 


9 = Wx A 

R + = R+ \ (V2) 


and extrapolate 7?~ from the interior. If outflow occurs at the boundary, there is only one 
ingoing characteristic (corresponding to /?+), and thus, set /?+ = extrapolate qp /?“, 

and s* from the interior. In the particular case of supersonic flow, all characteristics are ingoing 
if there is inflow, and are outgoing if there is outflow. Therefore, the dependent variables are 
specified with their free-stream values if inflow occurs, and extrapolation is used to determine 
the boundary flow variables if outflow occurs. 

At a distance far enough away from a 2-D lifting body, the lifting body can be viewed 
as a point vortex, with strength proportional to the circulation associated with the lift. The 
components of the induced velocity at the far-field boundary caused by the vortex can then be 
computed. Moreover, the effective velocity components at the far-field boundary are computed 
as (ref. 38) 


u — Ur K , cos q + / sin o 
v = nx si no — F cose; 


(5.3) 


where 
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n r ;) 

4t r R 
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Figure 5. Physical domain for airfoil calculations. 


Here, the subscript oc refers to free-stream values, a is the angle of attack, R and <p are the 
magnitude and angle of the position vector originating from a reference point at the body (i.e., 
quarter-chord point for airfoil) and extending to the far-field boundary point, respectively, c 
is the body length, and q is the lift coefficient. The polar angle <f> is defined as positive in 
the counterclockwise direction relative to a reference line (i.e., coinciding with chord for airfoil) 
emanating from the leading edge of the body and proceeding downstream. The Cartesian velocity 
components u and v of equations (5.3) are used to compute the local tangential and normal 
velocity components, respectively, required in the boundary conditions. 

Consider the case of a C-type mesh wrapped around an airfoil, and denote the outer boundary 
of a finite domain as Ei + (fig. 5). For airfoil computations, the boundary cells at Si are 
treated as described in this section. The boundary cells at S *2 are also treated in this way when 
the flow is in vise id. In the viscous flow problem, a portion of the boundary E 2 can generally be 
wake flow. If the boundary conditions applied at S 2 for in viscid flows are used for viscous flows, 
instabilities can occur. One way to treat the boundary cells at E 2 is to specify the pressure and 
extrapolate the variables /?, pu, and pv , which would satisfy the requirement of characteristic 
theory to specify one quantity. However, this approach results in pressure-wave reflections, which 
can seriously delay the convergence of the numerical scheme. 

An alternative boundary-point treatment is to extrapolate all dependent variables, allowing 
the outer and surface boundaries to determine a unique solution. Numerical experiments 
demonstrate that solutions obtained applying these two treatments are essentially the same 
near the airfoil. Furthermore, if the outer boundary is far enough away (i.e., 20 chords), there 
is generally no effect on global quantities such as lift and drag. The second approach shows 
noticeable improvement in the convergence behavior of the solution algorithm. 

For internal flows where the inlet Mach number is subsonic, the specified flow quantities of 
equations (5.2) are replaced with the total pressure, total enthalpy, and flow inclination angle. 
These conditions are usually known for internal flow problems. The Riemann variable R~~ 
is extrapolated from the interior of the domain. At a subcritica.1 exit boundary the pressure 
is specified, while the Riemann variable /?+, the total enthalpy, and the velocity component 
parallel to the boundary are extrapolated from the interior. If supersonic flow occurs at the 
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inlet or exit hound ary. the variables at that particular boundary are determined in the same 
manner described in this section for supersonic external flow problems. 

6. Basic Time-Stepping Schemes 

In section 6.1, the class of Runge-Kutta (R-K) schemes used for time integration is defined. 
The parameters associated with these R-K schemes and the requirements for determining the 
parameters are discussed. Then, stability analysis for the four-stage and five-stage schemes 
that are applied is conducted by considering a linear- wave equation. In section 6.2, stability 
properties of R-K schemes for systems of fluid dynamic equations are presented. This requires 
writing the Navier-Stokes equations in general curvilinear coordinates and defining associated 
Jacobian matrices. With this framework in place, an estimate for the time step is given in 
section 6.3. 

6.1. Runge-Kutta Schemes 

For problems where the area of a mesh cell is independent of time, the semidiscrete system 
of equation (3.2) becomes 

J t W <v + R(W;j) =0 (6.1.1) 

where R(W/j) is the residual function defined by 

R(w ij) = t~—~ (C c + C D + C AD ) W ( J (6.1.2) 


A variety of methods for t he integration of ordinary differential equations (ODE s) can be 
used to advance the solution of equation (6.1.1) in time. Single-step, multistage schemes 
(such a„s R-K schemes) are usually preferred, rather than linear multistep schemes (such as 
the Adams-Bashforth scheme), because multistep schemes require more storage and introduce 
implementation difficulties when combined with a multigrid method. A four-stage R-K scheme 
(ref. 1) that belongs to the class of standard RrK schemes and is fourth-order accurate in time 
Is used to solve a system of ODEs corresponding to the Euler equations. The four-stage R-K 
scheme can be written as 


W (0) = W ( n) 

w(E z= w ( °) — 

W (2) = W (0) - y-R (1) 

w (3) = w (0) - a/r (2) 


\ 


(6.1.3) 


W (4) _ w (0) - ~(R (0) + 2R (1) + 2R (2) + R. (3) ) 

_ ^y(4) 


where R^ = R.(W^)), the superscript n denotes the t ime level ?i A /. and the mesh indices 
(/, j) associated with the solution vector W are suppressed for convenience. If interest is only in 
steady-flow problems, then the higher order accuracy in time is not important, and other classes 
of multistage schemes can he considered. Schemes can he constructed with certain desirable 
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stability and damping characteristics that lead to efficient steady-state solvers. For example, 
the solution at the ( q + 1 ) t h stage (ref. 3) can be expressed as 


- w (0) _ 0?+1 a<R ( ^ 


where the residual function 


Riq) = q ( £ &V £ fW (r) + V £ D W (r) + ^DW< r > 


Vr=0 


r=0 


r— 0 


(6.1.4) 


(6.1.5) 


and the consistency conditions Ylfiqr ~ ^ ^2^qr — 1? and ^ 7 ?r = 1 must be satisfied. The 
basic parameters (where p = q + 1 (g = 0, . . . , rn — 1)) and the weighting factors /3q r , 
bqr , and yvyr must be prescribed to define the ??i-stage, time-stepping scheme. The desired 
stability and damping properties of the scheme provide the requirements for determining the 
basic parameters and weighting factors. Both hyperbolic and parabolic stability limits must be 
considered. The hyperbolic and parabolic limits are intervals along the imaginary and negative 
real axes, respectively, in the complex plane. The coefficients a p can be chosen to have the 
best possible hyperbolic or parabolic stability limit without special regard to the high-frequency 
damping characteristics of the scheme. However, if the scheme is used as a driver of a multigrid 
method, the scheme must effectively damp the highest frequency error components. 

Van der Houwen (ref. 39) gives the parameters a p that correspond to the maximum (or nearly 
so) attainable Cou rant -Fried rich s-Lewy (CFL) number. For schemes with an odd number of 
stages. Van der Houwen proved that the largest possible stability interval along the imaginary 
axis of the complex domain is (ni — 1). Vichnevetsky (ref. 40) conjectured that (m — 1) is also 
the optimal CFL number when w is even, and demonstrated this concept for 77? — 2 and 4. 
Sonne veld and Van Leer (ref. 41) proved that the (777 — 1) CFL number limit is valid when rn 
is even. Jameson (refe. 3 and 42) considers schemes with the a p s of Van der Houwen, and 
defines appropriate weighting factors for the artificial dissipation evaluations to yield a good 
parabolic limit. Although the a p s are obtained using only a hyperbolic stability limitation, 
they are still a good choice for a viscous flow solver, especially at high Reynolds numbers. That 
is, the convection (hyperbolic) limit on the time step remains the controlling stability factor for 
practical aerodynamic flows. 

Several members of the class of schemes defined by equations (6.1.4) and (6.1.5) have been 
analyzed in reference 42 by considering the model problem 


d w dw 9 d 4 w 

w + a— + , 4 Ax j-j 


= 0 


(6.1.6) 


Equation (6.1.6) is the 1-D, linear-wave equation with a constant-coefficient, third-order dissi- 
pation term. If the spatial derivatives in equation (6.1.6) are approximated with central differ- 
encing, then 


A/— = 1 - wj_ 1 ) - ( 4~( w j+2 ~ 4u’j+i + 6«>j* - 4w'j_ l + u>"_ 2 ) (6.1.7; 


where N — a At/ Ax is the Cou rant number. Taking the Fourier transform of equation (6.1.7), 
obtain 

dw 

At—r = ~ w tl (6.1.8) 


where the Fourier symbol 


dt 


z = —/iV sin 0 — 4f 4- — ( 1 — cos 0)^ 
a 


(6.1.9) 
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Here, i — y/— 1 , and 0 is the Fourier angle. If the residual function for any stage q is given by 

R {(l) = 0 C c + C ad )wM (6.1.10) 

then the amplification factor for an m - stage scheme is 

g(z)=l + f(9)z(0) (6.1.11) 

where 

f (0) — + K-2Z 2 ~h **■+ K m Z ni (6.1.12) 

Here, kj = a m with a m = 1 for consistency, and = ^/_i 1 with / = 2, 3... m. Since 
g(z) is analytic, the maximum modulus theorem guarantees that all contours \g(z)\ < 1 lie inside 
the absolute stability curve \g(z) \ = 1. For this subclass of schemes, which are schemes satisfying 
the requirement that \z\ < (rn — 1) (refs. 3 and 41), the optimal polynomials are defined as 


g(z) = P k (z) = i k ~ l T k _ 1 




(6.1.13) 


where Tf. is a Chebvshev polynomial, and k > 2. The coefficients K/ for the four-stage and 
five-stage schemes given in reference 39 are 


= 1 


k 2 - 


1 
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«3 = 


6 



an d 


= 1 


-2 = 2 


a:_‘> = 


3 

16 




i 

Kr> ~ 128 

respectively. 

In the more general situation, the amplification factor is not a polynomial in c . For example, 
consider the subclass of schemes defined by equation (6.1.5) that are called (???. n) schemes. The 
m refers to t he number of stages, and n designates the number of evaluations of the dissipative 
contribution. For example, assume a (m.2) scheme, where the numerical dissipation terms are 
evaluated on the first and second stages and frozen for the remaining stages (similar to the (4.2) 
scheme used). Let z r = 4i(~) and Zj — ?'T(r). Then, if w > 3. t lie / of equation (6.1.12) is 
replaced by 
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(6.1.14) 


f = K 1 - (Kiai^,. + K2Zi)z ( l + (K-2C*12r + «3‘i) 2 / 


+ ( — 1 ) >n (K, n -2 a lZr + K m — lZi)zj n 3 — (— l)' n K 


„ m — 2 


Consider the (5,3) scheme generally used, where the dissipation terms are evaluated on the 
first, third, and fifth stages, and frozen on the second and fourth stages. For this scheme, the 
weighting factors for the dissipation terms ( 7 ^. in eq. ( 6 . 1 . 5 )) are as follows: 


7 2o = r 3 
*'21 = 0 
7-22 = T 3 

7 30 = r 3 

731 = 0 , 

7 32 — T 3 
'33 = 0 
''40 = * 3*5 
''41 = 0 
^42 — 7 3 *5 

7 43 = 0 

7 44 “ 7 5 > 


(6.1.15) 


where I ;! = (1 — 7 ^), 1- = (1 — ~ . ) . 7 ^ = 0.56, and 7 = 0.44. The symbol of the time-stepping 
operator / for this scheme is given by 


/ = 


1 — a 4 *i( 1 


Q3-,:) - «4~3^l(«3~2~,- - 7.,-r) 


r 5''3-3~'- 


(6.1.16) 


where 

-1 = + “--/■ 

-2 = =i + T.j - r 

-3 = «aU - <n -i) 
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In a number of flow computations, four-stage and five-stage schemes are applied (refs. 19, 24, 
43, and 44). For these schemes, the residual function is 

R (,/) = f + £z>w( ()) +Xbv'- Cl/?W (r) j (6.1.17) 

By evaluating the physical diffusion terms on the first stage only, the computational effort of 
the scheme is reduced. This incompatibility with the computation of the numerical dissipation 
terms does not cause any deterioration in the performance of the schemes. The evaluation of the 
numerical dissipation terms on certain stages (and the weighting of these evaluations) provides 
two advantages. First, the parabolic stability limit can be extended, and the high-frequency 
damping can be improved. Second, the expense of calculating the dissipation terms can be 
reduced. Reference 3 provides additional discussion on the weighting of dissipation. 

The time-stepping parameters for the four-stage scheme are 

4 

1 

3 (6.1.18) 

1 

2 

04 = 1 > 

Since o ;/l _j = 1/2, the scheme is second-order accurate in time. The numerical dissipation terms 
are evaluated the same as for equation (6.1.14). For the model problem of equation (6.1.6), the 
absolute stability curve for this scheme is presented in figure 6(a). The hyperbolic stability 
limit (as determined by the extent of the stability interval on the imaginary axis) is The 

parabolic stability limit (as determined by the extent of the stability interval on the negative 
real axis) is 4. The dashed line (fig. 6(a)) represents the locus of the Fourier symbol as defined 
in equation (6.1.9), and must lie inside the \g\ = 1 curve for stability. Figure 6(b) shows the 
variation of the amplification factor (j with the Fourier angle 0. The scheme exhibits good 
high-frequency damping, which is crucial for a rapidly convergent multigrid method. When 
analyzing the stability and damping propert ies, it is important to include all components of the 
scheme. For example, if the stability limit of the algorithm is extended through the int roduction 
of an implicit residual smoothing procedure (discussed in section 7.2), some deteriorat ion in the 
high-frequency damping occurs (figs. 6(c) and 6(d)). 

In the case of the five-stage scheme, the basic parameters are 


<n 

«2 

°4 

nr, 



°1 = 
0'2 = 
= 
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(c) Stability curves with implicit residual smoothing; ( cl) Amplification factor with implicit residual smooth 

CFL = 4.8; ji = 0.6; k ( 2) = 0; = 1/32. ing; CFL = 4.8; ji = 0.6. 


Figure 6. Plots of four-stage R-K scheme; k ,2) = 0; k M| = 1/32; coellicients are 1/4, 1/3, 1/2, and 1. 





A very large parabolic stability limit is established by evaluating the artificial dissipation terms 
on the first, third, and fifth stages with the weightings of equations (6.1.15). Figure 7(a) shows 
that for this scheme, the hyperbolic stability limit is 4 and the parabolic stability limit is about 9. 
Figure 7(b) shows that this five-stage scheme also has good high-frequency damping. Figures 7(c) 
and 7(d) show the stability curves for this scheme when implicit residual smoothing is applied. 
In practice, with the implicit residual smoothing, a CFL number of 7.5 is used for this scheme. 
For the four-stage scheme, a CFL number of 5.0 is used. So, with one additional evaluation of 
the implicitly smoothed residual function, the CFL number increases by a factor of 1.5. 



(a) Stability curves with three evaluations of dissipa- 


tion: CFL = 3.0. 



(c) Stability curves with implicit residual smoothing; 
CFL = 6.0; d = 0.6. 



0 

(h) Amplificat ion factor with three evaluations of dissi- 
pation; CFL = 3.0. 



0 .5 1.0 1.5 2.0 2.5 3.0 3.5 


9 

(cl) Amplification factor will l implicit residual smooth 
mg; CFL = 6.0; d = 0.6. 


Figure 7. Plots of five-stage R-K scheme; — (J; k 4'U — 1/32: coefficients are 1/4. 1/6, 3/8. 1/2. and 1. 
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6.2. Stability of Rimge-Kutta Schemes for Systems 

The compressible Navier-Stokes equations can be classified as either hyperbolic-parabolic or 
incomplete-parabolic type (refe. 45-47). As discussed in section 6.1, these equations are solved 
numerically with a member of a class of multistage time-stepping schemes. By considering the 
hyperbolic (inviscid) and parabolic (viscous and numerical dissipation) operators separately, 
sufficient conditions for stability can be obtained. These conditions can be used as a starting 
point for estimating a time step for solving the full Navier-Stokes equations. 

To estimate the restriction on the time step due to convection and diffusion, we consider 
the two-dimensional, Navier-Stokes equations expressed in generalized coordinates (£,?;). These 
equations can be written as 


1) 

where W = J~ l W, with J 1 and W representing the inverse of the transformation Jacobian 
and the vector of conserved flow variables, respectively. The derivatives associated with 
the transformation from Cartesian coordinates (x,y) to (£,//) and the corresponding inverse 
transformation are related as 


dW <9F ( ^ <9f[,' ;) 

dt d£ drj dt drj 



li- 

1 

y>i 

-y* 

Ay 

ny. 

■ J- [ 




and the inverse of the transformation Jacobian is given by 

7 _l d(x,y) 

J =W^) =HV "^ r " V( = n 

The transformed inviscid flux vectors F«) and F<''> are given by 

F (?) = J _1 (^F ( ^ -^F (,/) ) 


( 6 . 2 . 2 ) 


F (,/) = J~ l (Vr F (0 - VyV { ' l] ) 


(6.2.3) 


Using the notation of reference 48, the viscous flux vectors f[P and f[P can be expressed as 




<9W 


, dw 

d£ 1 drj 

p('/) _ p(fj,<;) ^ W , 5(</.^) ^ w 

*’ drj d£ 


(6.2.4) 


(6.2.5) 


where 




a, d G {£•»/} 


( 6 . 2 . 6 ) 
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and //] ' ' is a viscous matrix obtained by rewriting the F r vectors in terms of primitive 
variables. The matrix M . which transforms nonconservative (primitive) variables to conservative 
variables, and the inverse of matrix M are given by 


M = 


-up 

-vp 


-1 


-1 


6* 


and 


A/ -1 = 


0 

P~ l 

0 

-(y-l)u 

1 0 


0 0 

0 0 

p~ l 0 

-(y-iyv ( 7 - 1 ). 


(6.2.7) 


(6.2.8) 


0 0 

u p 0 G 

v 0 p G 

H7-l ) -1 ^ 2 P u P v (7~ 1) _1 J 
respectively, with <t ft = (7 — 1 )(w^ + r^)/ 2 . The viscous matrices B^ Ji \ B^ J} \ and B^ 1 '^ 

are defined according to 


B ( a ' ij) = 


with 


IV 


B 1 '? = 

<1 h 


0 

0 

0 

W -2 


»-»0 , J 


B 


,a , 3 


Bj'i-U + B'y'jv + Bj-lyU B'^p - 1 


•CIO . d 
t , ‘ i . p(i, it 


0 

0 

0 


(6.2.9) 


yJ~)M f da d3 dad 3 


Rf 


x lh.db +l, dblh +flVa 


a,b 


( 6 . 2 . 10 ) 


B a ' d 


VjM ( 7 \ 

Re V7-1 ) 


—Vo • V/J 
Pr 


( 6 . 2 . 11 ) 


where the nondimensionalization of section 3.1 is employed, V is the gradient operator, and 
6 a L is the Kronecker delta, t’sing the Jacobian matrices .1 1 rl ) = dF^'^/d'W, where o G {«./?}, 
equation ( 6 . 2 . 1 ) can he rewritten as 


HI + th) “ Hi \ Hi + ihi j 

di) \ d£ dp 


( 6 . 2 . 12 ) 


where 


2 ('») - 


a y it - (7 - 1 )o yr 


0 o, 

-uUdd + o r 0 2 [/(") _ ( T _ 2)o,« 

— vU^'^ + a a <p’ atrv — (7 — l)o v i/ /d' 1 ) _ (- ; _ 2)a, l r 

irUdtf* -u,) a,* - (7 - 1 )«(/('*) aw - (7 - 1 


0 

(7-Do,- 
(7-1 )o.i/ 
7 / r(n) 


(6.2.13) 
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with U {a) — a x u + ctyV, and u t — 7 E — <f ? . Assume the coefficient matrices of equation (6.2.12) 
are locally constant in space and time, and transform to primitive variables. Then, obtain 


where 


aw P T ( (l »w P T (,„f)Wp 

~iT +A — + A — ~W 


+ (B 
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The eigenvalues of the matrices A * ^ are 

£/(«) 

£/(<») 

+ Cyjog + a] 

iW-CyJai+Ctl J 

The eigenvalues of t he matrices B ^ ” ° ' are 

0 

^ r > 

/iFi 

(2/1-f- A)rt J 


where 


r _ n/tM 1, 2 2 . 

r i — — 75 (°j- + «</) 

/re p 


(6.2.14) 

(6.2.15) 

(6.2.16) 

(6.2.17) 

(6.2.18) 
(6.2.19) 
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For the mat rices B 




, the eigenvalues are 

0 

v/7 M 


Re 

Re 


2-^Va V/?) 
rr p 


(U3/,) ,(Va-^)+M | Vo 


Re 


P 

(A + ‘.if i ) 

P 


(Vo • V0) - 


P 

(A + ft) 


|Vo| \Vi3\ 


( 6 . 2 . 20 ) 


Abarbanel and Gottlieb (ref. 49) showetl that the matrices of equation (6.2.14) can be simulta- 
neously symmetrized with the similarity transformation determined by 


0 0 0 

0 10 0 

0 0 1 0 

4 -pc 0 0 1 pc 

’ 4^ -l c -1 0 0 0 

0 10 0 

sr 1 = 

0 0 1 0 

rT^-rrr 0 0 

Using this similarity transformation, rewrite equation (6.2.14) as 


( 6 . 2 . 21 ) 


( 6 . 2 . 22 ) 


™ + A^— + A^— = + B^^r + (B^ + B^U) - 

cH di chf d (, 2 c)i ] 2 v / c 


,av 


,^v 


,5 2 V 


uY 2 V 


d' 2 V 




where 


A (f> ) = S~ l ~^ {a) S 
/?('*, d) _ s~ l ~B {a,,i) S 


(6.2.23) 


(6.2.24) 


with a, i) 6 {£,//}. 

Define as a discrete computational domain {(£ . ?/) : 1 < £ < L, 1 < 7/ < L, and Ad; = A i] — 1 } . 
where L = N, the number of mesh intervals in either coordinate direction. Let a discrete, 
vector function, such as V (/ Af . j A?/) = V (/._/) , be denoted by V;j. If the spatial derivatives 
of equation (6.2.23) are approximated with second-order cent ral differences, the semidiscrete 
representation 

j\ L j + C Vi J = 0 (6.2.25) 
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is obtained, where 

- n&nfaYij - (6.2.26) 

Then, taking the Fourier transform of equation (6.2.25) yields 

A^V/j = —AtCVij = ZVjj 
where the caret indicates a transformed quantity, and 

Z = Z C + Z D (6.2.27) 

with 

Z c =-iAt(A sin sin # /yl j 

Z^ — — At 4 sin" ) sin 6 ^ sin 0 rf +4 ,/,r ^sin 2 ^ 

2 2 

and ?* = \/— 1 . 

If all terms in the flow equations and the numerical dissipation are evaluated at each step 
of the R-K scheme, then the amplification matrix G for an 777-stage scheme is a function of one 
variable. In particular 

G(Z) — I T k\Z + k^Z^ + ■ ■ • + K ni Z ln (6.2.28) 

where k\ is a function of the coefficients of the R-K scheme. (See eq. (6.1.12) for definition.) Let 
A q{G) be any eigenvalue of G(Z), and let A q {C) be any eigenvalue of C. Also, let z — — A/A^(£). 
The eigenvalue A q{G) is related to z as 

A q{G) ~ 1 + T k^z^ + * * • + K m z m (6.2.29) 

The stability of a scheme requires that the amplification matrix satisfies the condition 
II G n || < C for all n. The spectral radius a of a matrix is defined as equal to the largest 
eigenvalue in absolute value. If the matrix G is normal (i.e., GG* — G*G), then its norm Is 
equivalent to the spectral radius cr(G). Thus, a normal matrix requires that <r(G n ) = a n (G) < C. 
This condition is equivalent to the Von Neumann condition for stability that requires 

a(G) < 1 (6.2.30) 

Hence, if Z and G are normal matrices, the condition for stability is 

|<7(~~)l = |s(-A<A,(£)) | < 1 (6.2.31) 

for all q. 

To determine sufficient conditions for stability, consider separately the hyperbolic and 
parabolic operators associated with equation (6.2.1). First examine the stability for the Euler 
equations (i.e., Re — oc in eq. (6.2.1)). The following theorem gives a sufficient condition for 
stability: 
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Theorem 6.1 Suppose the R-K scheme satisfies the property that 

[-* NcFL-i AV'F/J C {z e C : |.g(-)| < 1} 


(6.2.32) 


where Kppp is the CFL number for the R-K scheme. Assume smooth initial data , such that 
the Cauchy problem for the Euler equations is well-posed. If the Eule r equations are solved with 
this R-K scheme and second-order centered difference approximations for the spatial derivatives, 
the condition 


At [<r(a^ (€) + /M (,>) )] < N C fl 


(6.2.33) 


where |a| < 1 , and |/3| < 1 is sufficient for stability of the linearized problem. 

The proof follows directly. Since the matrices and ,4^ are symmetric, the amplification 
matrix G(Z) is normal. Furthermore, for a central difference scheme, the Fourier transform of 

the first derivative is a pure imaginary number. Hence, Z — —i At ^T^sin 0^ + sin 0 ^ , 
and so the result follows (remember that A£ = A?/ =1). 

Now consider the parabolic equation derived from equation (6.2.23) by eliminating the 
first-order spatial derivatives. The Fourier symbol of the difference operator is given by 
equation (6.2.27) when = 0. Let N j) denote the diffusion number of the R-K scheme (i.e., 
N d defines the stability interval along the negative real axis). Then, the maximum allowable 
time step is Np) A//), where A to = Aj^/(4 ft) for the 1 -D scalar diffusion equation with central 
differencing, with // being the diffusion coefficient. A sufficient condition for stability is defined 
as follows: 

Theorem 6.2 Suppose the R-K scheme satisfies the property that 


[wV D ,0]C {zeC:\g(z)\< 1} 


(6.2.34) 


where Np is the diffusion number for the R-K scheme. Assume smooth initial data , such that 
the Cauchy problem for the viscous equations is well- posed. If the viscous equations are solved 
with this R-K scheme and second-order centered spatial differencing, the condition 

max (At a p) < Np ( 6 . 2 .35 ) 


with- 


er p —4a 


sin 2 - 7 -+ s \ n 2 JL L^giCn) _j_ ) s j n Qe sin 0 tj 


is sufficient for stability of the linearized equation. 
Since 


Z = -At 


4 



^ J ^sin 2 -jj-+ fii'/'O) sin^ 



sin Be sin 0 {} 


(6.2.36) 


is a. negative real symmetric matrix, G(Z) is again normal. Therefore, the proof is similar to 
the proof of theorem 6.1. Note that if the cross-derivative terms are neglected, the inequality of 
equation (6.2.35) reduces to 


At 


4<t(b {U) + 


< N 


D 


(6.2.37) 
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Next, suppose numerical dissipation in the form of fourth-difference terms is added to the 
discrete viscous equations. Then 


£V/j = 

+ e l4) <r(AW)6£Vij + £ i4) a(A {ll] )6 4 V ;j 


(6.2.38) 


and the Fourier symbol of equation (6.2.27) is replaced by Z — Z]j + %AD' where Zq is defined 
the same as in equation (6.2.27), and 


Zad =-16AF- (4) a(.4 


i'«sin 4 |- + 


it 7 ?) ■ ^ 
v " sin 


■i) 


with / being the identity matrix. 

Lemma 6.3 If only the viscous terms and the fourth-difference dissipation terms are considered , 
then the condition 

max [At (a D + <Jad)] < N D (6.2.39) 

with <J [) the same as given in equation (6.2.35) and 

a AD = IG^a^si^f+A^hm 4 

is sufficient for the linearized stability of the multistage scheme. 

The proof is the same as the proof for theorem 6.2. If cross-derivative terms are neglected, the 
inequality of equation (6.2.39) becomes 


At 


4cr(£(U) + + < N d 


((5.2.40) 


6.3. Time Step Estimate 

Now consider the situation where the R-K scheme simultaneously satisfies the properties of 
equations (6.2.32) and (6.2.34). Also consider the 1-D case in the f direction. The R-K scheme 

then depends on the matrix Z , where Z = — At sin 6 + 45^^)sin 2 0^/2^ . Since A^ 

and do not commute, the matrix Z , and thus the amplification matrix G\ are no longer 

normal matrices. Hence, the Von Neumann condition on the largest eigenvalue of matrix G is 
now a necessary, but not a sufficient, condition for stability. Thus, there is no simple way to go 
from properties of matrix Z to properties of matrix G. The spectral mapping theorem relates 
the eigenvalues of matrix Z to the eigenvalues of matrix G. Since the eigenvalues do not tell 
the entire stability story, energy estimates based on norms must be used. However, no simple 
relationship exists between the norm of matrix Z and the norm of matrix G. 


In practice, a simplified stability condition is used to estimate the time step. There is no 
strict mathematical proof of stability with this condition; nevertheless, it seems to work well. 


Oonsi der 


At C = 


n cfl 

<?C 


At 


D 
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where 


<JC = <t(A^) + <r(A^ ) 

(tjj— <t(B^£^) + a(B^ rjJ ^)- h <t (S^ T ^) 


1 1 

A/c* + A//> 


so that 


ff C- + ( N CFL /N D ) ° D 


(6.3.1) 


Schemes in which NcFL ~ ^ r D have been considered. The time step for each cell in the 
computational domain is then computed as 


where Q is the cell area and 


NcfJ. 
Ar + \d 


\c — Ac -f A ii 

= (^d)( + (^d)i; + (^d)(i 


(6.3.2) 


with A^ and X f/ defined by equation (4.2.5), and 

\ \fi M If 1 ^- 1 / 2 , 2 , 

(Xd) ^-rTJp^ q {r " + y » ] 


{Xd) '^rT7p 7 n_1( *l + ^ ) 


(6.3.3) 


—(A + + (A + t l )\J{yq + x ) ;)(y| + ^ ] j 

For the thin-layer, Navier-Stokes equations, take Xp = (A^),,, where ?/ is the direction normal 
to the boundary layer. 

Remark 6.4 So far , only central differencing is considered for the spatial approximations in 
estimating a time step for an explicit R-K scheme. Since a numerical-flux function for an 
upwind scheme can generally be expressed as the sum of a centered (physical) contribution and 
a numerical dissipation contribution , then equation (6.S.1) is also a reasonable estimate for the 
time step when an upwind scheme is used. 

7. Convergence Acceleration Techniques 

7.1. Local Time Stepping 

The first technique employed to accelerate convergence of the basic explicit time-stepping 
scheme to a steady-state solution is local time stepping, where each cell is updated using an 
individual time step. For simplicity, the one-dimensional Euler equations are used to understand 
the meaning of local time stepping from the discrete point of view. Suppose the Euler equations 
are written in the form 

r)W a W _ 

dt + dx 

If the equations are discretized in an explicit sense, then 

AW = 4 W" 
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where the tilde indicates that the vector or matrix is for the complete discrete system. The 
block matrix ,4 is a function of At n ^ n /Ax, with A t mm being the minimum time step permitted 
in the domain, and each element of ,4 being a 3 by 3 matrix. Let the explicit matrix >1 for the 
system of difference equations be preconditioned by the diagonal matrix A, given by 

A = Diag [ A t^I A/3/ ••• At N _. 2 I A t N -\l] 

where 

: = -Ai_ i = 2 N - 1 

A^min 

Here, A/,; is the largest local time step allowed by stability, and I is a 3 by 3 identity matrix. 
This process results in a significant speedup in the t ransport of information, and an increase by 
roughly a factor of two in the convergence rate of explicit schemes. 

7.2. Residual Smoothing 

The local stability range of the basic time-stepping scheme can be extended by applying 
a procedure called implicit residual smoothing. This technique was first introduced by Lerat 
(ref. 50) for the L ax- Wen droff scheme, and later devised by Jameson (ref. 51) for R-K schemes. 
The constant-coefficient approach of Jameson is discussed in section 7.2.1. Some basic properties 
of residua] smoothing are also presented. Then, variable coefficients for implicit smoothing are 
discussed. The coefficients introduced in this paper, and those of Martinelli (ref. 27), are derived 
and compared. In section 7.2.3, coefficients are developed for implicit residual smoothing that, 
allow a time-step estimate independent of a physical diffusion limit. 


7 .2.1. Constant coefficients . The constant-coefficient, implicit residual averaging of 
Jameson (ref. 51) can be applied in two dimensions using the factored form 


( 1 - ^ Ac )(1 - /?,, V,, 


(7.2.1) 


where the quantity V A is a standard second-difference operator, and thus 

T-7 a ^r( "0 _ TT ( »«) rwT ( ”0 , TT ( '" ) 

The quantity l3 is a smoothing coefficient, and (£,77) are the coordinates of a uniformly spaced, 
computational domain. The residual of the unsmoothed scheme is defined bv 


'rd m ) _ n d r r \AT^ nt ^ ) 1 r 

■ij - a ' n ft. 


i[£ c .w'"( 1 ] + C D W?J + AD ( "'>] (rn =1,5) 


(7.2.2) 


and computed in the Runge-Kutta stage ?n, AD( m > is the total artificial dissipation at stage 
777, and is the final residual at stage 777 after the sequence of smoothings in the £ and 7/ 

directions. A tridiagonal system of equations is solved for each coordinate direction to obtain 
the unknown residuals 7 To determine /?£ and f3 t ^ Jameson (ref. 51) considers the model 
problem of equation (6.1.6) without numerical dissipation (i.e., the convect ion equation). Then, 
the semidiscrete equation (6.1.7) becomes 


*£= -A iCwj 

at 
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with 


At Cu!j = y ( Wj + j - Wj_ , ) 

and the Fourier symbol of the difference operator — AtC is given by 

z — — AtC = —i N si n 0 


(7.2.3) 


Let AT Fj define the correction, or residual, obtained from the implicit smoothing procedure, so 
that 

(1 -/JVA) A wj = A wj (7.2.4) 

and the Fourier symbol of equation (7.2.3) is replaced by 


A' sin# 

1 + 4/3 sin" % 


(7.2.5) 


A sufficient condition for stability is as follows: 


max b| < jV 


7.2.6) 


for all 0, and N* is the Courant number of the unsmoothed scheme. Solving for sin 0 and cos 0 
corresponding to the maximum of |^ | yields 


sin 0 = 


/I + 4/3 
1 + 2/3 


1 + 2/3 


(7.2.7) 


Using equations (7.2.7) and the sufficient condition of equation (7.2.6), the smoothing coefficient 
is determined by 


1 ( N 


^4 - 


1 / At 


4 A t 


(7.2.8a) 


(7.2.8b) 


where At* is the time step of the unsmoothed scheme. In subsequent discussion, this i3 will be 
referred to as the 1-D smoothing coefficient and will be designated by ;di-|) 

Instead, consider the diffusion equation 

d w d 2 w 

~di = 

If the spatial derivative is again approximated with a central difference, and a Fourier transform 
Ls taken of the resulting semidiscrete equation, the Fourier symbol of the product of A / and the 
difference operator are given by 

- = -Arsin ' 2 

where the diffusion number A'/j> = 4 Ai/ij Ax^ . If the residual smoothing operator of equa- 
tion (7.2.4) is applied, then 

, = (7.2JI) 

1 + 4/?yy sill ^0/2 
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A sufficient stability condition for an R-K scheme is 


max |~| < Nq 

for all 9 , where N* D is the diffusion number of the unsmoothed scheme. With the same procedure 
employed for the convection equation, the smoothing coefficient is determined as 



Thus, for the scalar diffusion equation, the smoothing coefficient 0jy is proportional to 
fiAtf(Ax)^. As will be shown in section 7.2.3, this type of 0 can be combined with the type of 
/ 3 given by equations (7.2.8) to yield a formulation suitable for a convection -diffusion equation. 

Some properties of implicit residual smoothing are now examined. If residual smoothing is 
applied on each stage of an R-K scheme, the stability function given in equation (6.1.11) still 
applies, with the ^ of the original (basic) scheme modified as in equation (7.2.5) or (7.2.9). This 
stability behavior leads to the following theorem: 

Theorem 7.1 Let C be the Fourier symbol of any discrete spatial operator for the convection- 
diffusion equation. Let equation (6.1.11) be the stability polynomial for an explicit m- stage R-K 
scheme . Apply implicit residual smoothing, as in equation (7.2.4), after every stage of the R-K 
scheme . 

If the original scheme is unconditionally unstable , then the smoothed scheme is also uncon- 
ditionally unstable. If the original scheme is conditionally stable, then the smoothed scheme can 
be made unconditionally stable by choosing 0 At sufficiently large. 

Proof. Let 2 and z s be the symbols of the original and smoothed schemes, respectively. Then, 


1 + 4/? sin 2 0/2 


(7.2.11) 


Define r as the position vector corresponding to r. Thus, r emanates from the origin of the 
complex domain and has magnitude \z\. Let r b be the position vector associated with z#. If 
the original R-K scheme is stable, then r does not terminate outside the stability region 5, 
determined by the Von Neumann condition \g(z)\ < 1, where g(z) is the stability polynomial. 
If the original scheme is uastable, then there is no At > 0 small enough to allow r to be in 
.S'. Since the denominator of equation (7.2.11) merely acts as a scaling factor of r, the residual 
smoothing cannot stabilize the unstable original scheme. 

Suppose 3 is proportional to Ai, as in equations (7.2.8) or (7.2.10). Then, r s does not 
terminate outside the boundary of S for any value of A/, since 0 can always be made sufficiently 
large. Moreover, the scheme is unconditionally stable. 

Remark 7.2 Even though the explicit R-K scheme can be made unconditionally stable with the 
implicit residual smoothing, th ere is a practical limit on the time step when solving the hyperbolic 
problem and taking 0 oc At 2 , as in equations (7.2.8). That is, if At is too large, convergence 
slows down. 

Lemma 7.3 Apply an explicit m-stage R-K sch eme with implicit residual smoothing to the scalar 
equation 


where £* 4) is a constant coefficient . Assume that the residual smoothing coefficient j3 is 
proportional to N 2 . Then, the symbol of the smoothed scheme vanishes as At — oc and the 
stability polynomial g(z ) — 1. 

Proof. The symbol for the difference operator of equation (7.2.12) when implicit residual 
smoothing is applied is given by 


— Ifie^jVsin 4 0/2 
1 +4(3 sin 2 0/2 


Using 0 ~ jV 2 and taking At to be large, then 


4~ (4W 


0 

2 


Therefore, z# — * 0 as At — oc. From equation (6.1.11). it follows immediately that </(;:) — 1 as 

c,-0. 

Remark As evident from Lemma 7.3 , the limit on the extetision of stability with the implicit 
smoothing arid equations (7/2.8) is caused by the requirement to have a certain background 
dissipation (i.e.. high-frequency damping). If 0 ~ N . as in the parabolic problem, then the 
symbol z y does not vanish. 

The use of constant coefficients in the implicit treatment (eqs. (7.2.8)) proves satisfactory 
(extending the Oourant number by a factor of two to three) even for highly stretched meshes 
of viscous-flow computations (ref, 16), provided additional support such as enthalpy damping 
(ref. 1) is introduced. However, the use of enthalpy damping, which assumes constant total 
enthalpy throughout t he flow field, precludes the solution of problems with heat-transfer effects. 
By using variable coefficients 0c and 0 ff , which account for the variation in mesh -cell- aspect 
ratio, residual smoothing can be applied without the support of enthalpy damping. 


7/2.2. Variable coefficients. The alternating direction implicit (ADI) scheme and the im- 
plicit scheme of Lerat (ref. 52) exhibit a functional dependence of variable smoothing coefficients 
on the characteristic speeds Ac and A,,. as defined in section 4.2 of this report . Appendix B shows 
this functional dependence. Then, with a 2-D stability analysis similar to the IT) analysis of 
the previous section, variable smoothing coefficients can be obtained as 


0 c — max 
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1 + r nt 
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(7.2.13) 


where again r f/ c = X^/Xci. The limiting cases are 0c — th-Th ^ as 7 7/( — ’ 0 and 

0c — 0, 0 tj — 0 1 — 13 . as r (j c — oc. 

A problem exists with the smoothing coefficients of equations (7.2.13). In the typical case 
of N/N* — 2. the smoothing coefficients vanish when r /f c — 1, making the scheme unstable. 
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Martinelli (ref. 27) eliminates this difficulty by modifying the residual smoothing coefficients of 
equations (7.2.13) as follows: 


where 


0 c = max 


0, f = max 


4 

l 

4 



<t>( r ,£ ) 

1 + r 'i£. 



N 




N* 1 + r 


-1 



,0 


> 


(7.2.14) 


<p(r) = 1 + r** 


(7.2.15) 


and the exponent ( is taken to be 2/3. The function <p Ls the same function used for scaling the 
artificial dissipation coefficients. The introduction of this function seems appropriate because 
of the direct relationship between residual smoothing and artificial dissipation. For example, a 
desired high-frequency, damping behavior of the scheme can be maintained when the dissipation 
is increased by increasing the residual smoothing. 

The variable smoothing coefficients 0 £ and 0 , } of equations (7.2.14) cannot be uniquely 
determined from a sufficient condition for stability, as the constant coefficient 0 was in 
equations (7.2.8). Wigton and Swanson (ref. 53) use an additional constraint to derive the 
parameters of equations (7.2.14). For completeness the short derivation of reference 53 is 
presented. 


Consider the following sufficient condition for stability: 


1 JV 1 1 

N* 1 + v /TT4/T + N* 1 + r~} ^/TTW, ~ 1 


(7.2.16) 


as derived in appendix B. Let. the Courant numbers for the two coordinate directions £ and // 
be given by 
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A A<c + A tf 


A<Hct = Af A,, 

Airy A^ + A II 


N ' 
1 + r >£ 
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1 + r 


-1 
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(7.2.17) 


where A/ ac t is the 2-D allowable time step for convection, and At c and A/ r; are the corresponding 
1-D time steps. If the Courant number A r in equation (7.2.16) is replaced according to 
equations (7.2.17), the result is 


A'£ 1 

N* y/Y+Wt 


1 

x* yr + ah,, 


< i 


(7.2.18) 


As suggested by lemma 7.3, 0 ^ and 0 , } should be as small as possible and still maintain stability. 
With this objective in mind, 0 c + 0 n is minimized subject to equality in equation (7.2.18). A pply 
the method of Lagrange multipliers and consider the function 


F{ i A/) — 0c + 0i) + & 


( N t 1 . Ay 1 

V V * s/TT4JJ7 N* yi+4^ 
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After equating the partial derivatives of F ( ,+ . ) . with respect to lie and 3,j. to zero, obtain 


1 + 4 //^ ) 


3/2 


= N 


T ' ( 1 + Wnf /2 


(7.2.19) 


Now, solving for ^/1 + 4 i3 t) and substituting the resulting expression into equation (7.2. 16) yields 


JL 1 1 [, , (XY\ 

N* 1 + yji + 4/+ ^ 1 + r~} V y V J 


1 / 3 - 

< 1 


(7.2.20) 


With the equality of equation (7.2.20) holding, and using equations (7.2.17). obtain the 
smoothing coefficients of equations (7.2.14). As shown in reference 53, the function O arises 
without any consideration of numerical dissipation terms. The role of the o function is to 
connect the values of /? corresponding to low-aspect-ratio and high-aspect-ra tio cells. 


The variation of ji^ from equations (7.2.14) with r 1f c is shown as a solid line in figure 8(a). 
For this curve, the ratio of Gourant numbers N/N* is assumed to be 2. Observe that lie — * 3\-D 
for j‘ ff £ — 0, and fig — 0 for r f g — * oc. In the case of i* lf e = 1, 3^ = 3 r f — ( a va l ue °f 0.75 

when N/N* = 2). Based upon numerical calculat ions for inviscid flows using t ypical in viscid 
meshes, smoothing coefficients 3c and 3ij that are constant with a value of about 0.4 result in 
rapid convergence. Values for lie > 0.75 and/or 3 f f > 0.75 can cause a. significant slowdown 
in convergence. To provide improved smoothing coefficients, when % 1 the formulas of 
equations (7.2.14) can be replaced with 
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(7.2.21) 


where V’ is a parameter to be specified. Here, the connection function O is removed by introducing 
V’. Figure 8(a) shows the curve representing //c when N/N* - 2 and + = 0.25. For the case of 
r, K c = 1 , /+ is 0.39. 
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Figure 8. Implicit residual smoothing coellicient lor two forms of variable coefficients. 
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To investigate stability using the smoothing coefficients of equations (7.2.21), consider the 
sufficient stability condition 


N 1 sin (h N 1 sin 0„ 

1 -j 1 < i 

N* 1 + r,£ X£\ n N * 1 -f- r~£ \(\>i 
for all and 6 rf . In equation (7.2.22) 


= 1 + 2/% (1 - cos ) 
\ n = 1 + 2j3,i ( 1 - cos0, ; ) 


(7.2.22) 


(7.2.23) 


and = \ rj /\c. This condition comes from the 2-D stability analysis given in appendix B. 
(See eq. (B18).) If r 7 ^ <C 1, the condition of equation (7.2.22) reduces to approximately 


N 

N* 


sin 0c 


< i 


for all (F . Using equations (7.2.7), obtain 
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N* 
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x/ 1 + 


< 1 


(7.2.24) 


(7.2.25) 


Substituting for /% according to equations (7.2.21 ). the inequality equation (7.2.25) is satisfied. 
Now, if r (/ c 1. the condition of equation (7.2.22) reduces to approximately 


N sin On 

~ \i) 


< 1 


(7.2.26) 


for all 0 t) . From equations (7.2.7) and the definition of /% in equations (7.2.21), the inequality of 
equation (7.2.26) is satisfied. Consider the case of r^c — 1. Assume Oc = 0, t and Jc = 3 tj . Then 
equation (7.2.22) becomes 

F < 1 (7.2.27) 


for all tic, where 

^ iV si n Oc 

F= ~ — 

It can be shown that 
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Then 


N 1 

Fn,ax - °' 9 ~ yi + 4 /% 


(7.2.28) 


By substituting for 3c and taking t’ ~ 0.11, the condition of equation (7.2.27) is satisfied. From 
numerical experiments, this estimate for V’ seems to be conservative. A value for of 0.125 
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worked well for both central and upwind schemes. Moreover, for central schemes V’ = 0.25 also 
works without causing stability problems. 

Calculations were performed for a case of inviscid, transonic flow 7 over an RAE 2822 airfoil to 
evaluate the variable smoothing coefficients of equations (7.2.14) and (7.2.21). A typical inviscid 
mesh with 224 by 32 cells w r as used. Figure 8(b) shows the convergence histories corresponding 
to the formulations of Martinelli (eqs. (7.2.14)) and Swanson (eqs. (7.2.21)). Convergence is 
measured by the logarithm of the rootr me an -square of the residual of the continuity equation. 
For each computation, the basic explicit scheme of equation (6.1.4) and a multigrid method 
(described in section 7.3) were employed. The average rate of reduction of the residual was 
defined by R j — (rate)A where R is the residual for the continuity equation, the subscripts / 
and i mean final and initial values, respectively, and N denotes the number of multigrid cycles. 
With the coefficients of equations (7.2.14), the average rate of residual reduction is 0.889, while 
with the coefficients of equations (7.2.21), the average rate of residual reduction is 0.789. As 
expected, the two formulations exhibit only small differences in convergence behavior in the case 
of turbulent flow calculations, since the high-aspect-ratio cells of the mesh, usually defined to 
resolve the boundary layer, determine the convergence rate. 


7 .2.3. Removal of diffusion limit . The diffusion restriction on the time step (eq. (6.3.1)) 
can be a significant factor in viscous regions of a flow 7 field, causing excessive restrictions 
on the allowable time step At. In this section the diffusion-based, smoothing coefficient of 
equation (7.2.10) is utilized to construct a new smoothing parameter that allows the removal of 
this diffusion restriction. 

Considering the thin-layer form of the 2-D Navier-Stokes equations allows use of the 
smoothing coefficient of equations (7.2.21) in the s t ream wise- 1 ike (£) direction. A possible 
formulation for the normal tj direction depends on a diffusion- type f3 near the surface, and 
a convection-type /? when the viscous effects are no longer important. Consider the dependency 
in equation (7.2.10) on the ratio of the actual At to the At of the basic explicit scheme. To 
remove the diffusion limit on the time step, the actual time step must be independent of diffusion 
effects. Thus, set N jj = 0 in equation (6.3.2), giving 


At = At act 


NQ 

A^r + Xfj 


(7.2.29) 


where Q is the area of the mesh cell being considered. In the part of the boundary layer wdiere 
diffusion effects dominate, define the time step of the unsmoothed scheme ( A/*) by 


where 


At* = (A t D ) fl 
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(7.2.30) 
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If (A i o)if — A / a ct i &d ~ 0- This means that, the full parabolic stability limit is being used for 
the physical diffusion terms. Since the numerical dissipation of the scheme is not included in 
the analysis, replace equation (7.2.31) with 


(0d)i) 


\_ 
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~^ar I 

(A^)„ 



(7.2.32) 


where C\ is a constant. Equation (7.2.32) accounts for any possible influence on the stability 
caused by a single evaluation of the physical viscous terms in the multistage time-stepping 
scheme. By using equations (7.2.29) and (7.2.30) to replace A/ ac t and (A /£>),;, respectively, 
equation (7.2.32) can be rewritten as 




C\ 


N (A p) rf 

I\i D 



(7.2.33) 


which can be approximated by 


(/?£>)// — 


J 7-7 (A p)i) 
4 1 Ac T A , f 


(7.2.34) 


Either equation (7.2.33) or equation (7.2.34) can be used, provided the constant is defined 
properly. Both equations successfully remove the diffusion restriction on the time step. In this 
paper, the simpler form of equation (7.2.34) is used, and has also been considered by_Radespiel 
and Ivroll (ref. 54). Numerical experiments have shown that a satisfactory value for C\ is 5. 


For the full Navier-Stokes equations (including all viscous terms), a coefficient ftp for the 
streamwise-like direction (£) should also be defined. Using the form of (/?/})// hi equation (7.2.34), 
(&d)( * s defined as 
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where 


(*ok = 


Vr M 2H. 

Re pPr 


l ( x i + y'n) 


The variable coefficient of equation (7.2.34) generally cannot be used alone. For example, in 
an airfoil flow, goes to 0 too fast at the leading edge, resulting in a 0-value in the in vise id 

region. This difficulty is overcome by calculating as 

0i) = max[( /?£>), 

where is defined by equations (7.2.21). In a similar manner, (/?£>)£ in equation (7.2.35) is 

redefi ned. 

According to the theory presented in reference 55, the residual smoothing is evaluated only 
on the even steps of an H-K time scheme. In practice, the residual smoothing is evaluated during 
every stage, which is more expensive but produces a more robust algorithm. 


7.3. Multigrid Method 

The concept of multigrid acceleration of an iterative scheme was first suggested by Fedorenko 
(ref. 56). The fundamental ideas of this approach currently used in many applications are 
principally due to reference 57. Although most of the theory developed for the multigrid method 
is for elliptic problems, a number of effective multigrid solvers (refs. 3, 58. 59, and 60) have been 
constructed for the Euler equations of gas dynamics, which are hyperbolic. Transonic and 



subsonic flows have been computed with these solvers. Some multigrid methods (refs. 19, 20, 
22, and 23) have also been devised for the numerical solution of the compressible Navier-Stokes 
equations. In section 7.3.1, the basic theory of the multigrid process is briefly reviewed. Then, 
the operators used in the present method are defined. Section 7.3 concludes with a discussion 
of the various elements of the multigrid technique of this work. 


7 .3.1. Basic concepts of multigrid methods. In the simplest sense, the multigrid method 
involves applying a sequence of grids to solve a discrete problem. More specifically, a faster 
rate of development, of the solution on a fine grid is achieved by approximating the fine-grid 
problem on successively coarser grids in the sequence. With suitable coarse-grid approximations 
of the fine-grid problem, the low-frequency error components on the fine grid appear as high- 
frequency error components on the coarser grids. The low-frequency components on the fine grid 
where the discrete solution is desired are precisely the error components that dramatically slow 
the convergence of single-grid schemes. Thus, with a good high-frequency damping scheme, an 
effective multigrid process (i.e., much more rapid removal of low-frequency errors than a single- 
grid scheme) can be constructed. As will become evident, the driving scheme for the multigrid 
process is not only important for providing smoothing on each grid, but also for removing high- 
frequency errors resulting from interpolation of corrections for the fine-grid approximation. 

Two additional advantages are derived from displacing part of the effort in solving a set of 
discrete equations to coarse grids. One advantage is that the larger mesh spacing permits larger 
time steps, meaning that information is propagated rapidly in the domain of interest. Moreover, 
for explicit time-stepping schemes such as the multistage schemes described previously in this 
report, the increased time step Is particularly important because the allowable time step depends 
on the speed of sound. This acoustic dependence is even more critical for viscous flows. A 
second benefit of the coarse grids is that they require less computational work. For example, 
in two dimensions, the computational effort needed is decreased roughly by a factor of four on 
successively coarser meshes. Thus, the objective of the multigrid process is to spend much more 
time on the coarse grids than on the fine grid. 

The basic ideas of the multigrid process are revealed by considering the continuum problem 

CW(x , y) = S(x,y) 

AW (a, y) - ${x, y) 


where the first equation Is associated with the domain Q, and the second equation is associated 
with its boundary 7 dQ . The symbols C and A are general nonlinear, differential operators, and 
both S and <F are source terms. Let Go,Gi, ...,Gy be a set of grids, where Gy is the finest 
grid, and each successively coarser grid Gfc(k < N — \ ) is generated by eliminating every other 
mesh line in each coordinate direction of the next-finer mesh. The discrete problem on Gy is as 
follows: 

AvWy(ar,y) - S N (x,y) (x,y £ Gy) 1 

Ay Wjy (x,y ) = T N ( x,y £ dG N ) J ) 

and IV y is the exact discrete solution. If w\<(x,y) is an approximate discrete solution, 
equations (7.3.1) can be written as 


Cswj v = ,$'y + Rj y 
Ay U y 3 = <J>y + (Rq)n 


(7.3.2) 
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where R y and (Rb)n are residual functions. Subtracting equations (7.3.2) from equa- 
tions (7.3.1) gives the residual equations for the G T y problem. That is, 


C \ liy — C \ riv \7 — —R \7 

Ayr W\r — Aywy =~(R b ) n 

These equations can be adequately approximated on Gy_i if the residual functions and 
corrections ( lUy — u: y) are smooth. Smoothing is accomplished by performing an iteration 
with an effective high-frequency damping scheme. The approximations of the residual equations 
on the coarser grid G\_ \ are 


£a _i W'v-i - £(lg~ l w N ) = -I%~ l Rx 
A.v- , W', v _ ! - A (I*~ l w N ) = -I*~\r b ) n 


(7.3.3) 


where 7y 1 Ls a restriction operator. Note that if R,\ = 0, then IVv - 1 = /y 1 w ,v . and once a 
steady state is reached on the fine grid, all corrections on the coarse grid are 0. Furthermore, for a 
linear problem, the two terms on the left-hand side of equations (7.3.3) can be combined and the 
error equation £(error) = — /y -1 /£yr is obtained. In general, the operator Ij n is used to indicate 
restriction when / > m and prolongation when / < m. Thus, a restriction operator transfers 
information from a fine grid to a coarse grid, and a prolongation operator (i.e., interpolating 
polynomial) transfers information from a coarse grid to a fine grid. Equations (7.3.3) can be 

rewritten as 

£,V_1 W ’.V-1 = %- I 


where 


A.v_i«’.v_i = 

S’.v-i = R n~\ + Rn- i 


(^b).V- 1 - + ( F b) ,V_1 
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r 1 T A — 1 

ejV-1 = 1 N 


(Fb)n- 




{—R.x ) + Cv-i(/.v Vv) 


( 7.3.4) 


Thus, the discrete problem on Gyr _ \ has the same form as that on G \ , except the forcing 
functions of equations (7.3.4) are added to the residual functions. An improvement to the 
approximate solution w\; can be obtained by adding a coarse-grid correction. The fine-grid 
solution is then given by 


W’iV + - tv '"’.V 


A'-l. 


where the correction (re.v— l ~ ' / y * w .\ ) is an approximat ion to the smoothed function M y — tr y . 

In this work, the smoother chosen to solve equations (7.3.3) is a multistage R-K scheme of 
the type discussed in section (5. Thus, a time derivative is added to the steady-state equations, 
and t he resulting equations are advanced in pseudotime with several iterations of the multistage 
method. Usually, one complete R-K time step is performed on the finest mesh, and two or three 
time steps are performed on coarser meshes. 

Instead of immediately passing a correction from G\_\ to G y. the solution u y-i and 
residual /f\_i can be restricted to the grid (7\-*2- Iteration sweejxs can then be performed to 
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obtain a smooth approximation of the correction function If ^v-2 ~ w N-2- If this correction is 
passed to i, iterations are performed, and the correction of G\t_± is transferred to Gy , a. 
multigrid cycle of three grids is completed. This cycle is called a V cycle. 

There are other fixed-cycle strategies (i.e., W cycle), and variable- cycle strategies that 
depend on a prescribed residual level, or a certain slowdown, in smoothing rate before changing 
to a coarser grid problem (ref. 57). For each coarse grid Gj~ in a cycle, the full current 
approximation w k and the initial (basic) approximation wk + 1 (the approximation on grid Gfc+i 
that is transferred to grid G j.) are stored. The approximation w ^ is the sum of the G correction 
and the basic approximation. Brandt (ref. 57) refers to a scheme that stores the full, current 
approximation rather than only the correction as the full approximation storage (FAS) scheme. 


7.3.2. Ti^ansfer operators. The intergrid transfer operators employed in the present 
multigrid method were introduced by Jameson (ref. 3) and assume that the unknowns are stored 
at the center of a mesh cell. The restriction operator for the residual is defined by 


Jfc+l^fr+1 “ (^Ar+l^+l)/ (7.3.5) 

tk /=! 

where the residual function Rfc+i is expressed in the usual way as 

R k+\ ~ o A*+1 u T+1 

with and + i denoting the spatial- discretization operator and the cell area, respectively, 

on grid G^+j. Thus, the modified residuals of the four fine-grid cells corresponding 

to a coarse-grid cell are summed. In this manner, the residual transfer operation is conservative. 
To transfer the solution from G to G*., the following volume-weighted operator is used: 

y^4 

lk+im+i= ^^4 

Again, the summations are over the four fine-grid cells, and the operator conserves mass, 
momentum, and energy. The prolongation of corrections from G/. to G^+i is accomplished 
with bilinear interpolation. 

In elliptic multigrid methods, the residual- restriction operator is frequently defined as the 
adjoint of the correction-prolongation operator, meaning that one operator is the transpose of the 
other. (See appendix C for discussion of the adjoint property.) Such a relationship is convenient 
for analyzing multigrid schemes (ref. 61). In typical multigrid methods using a cell-vertex, finite- 
volume formulation for spatial discretization (refs. 3 and 24), the restriction operator is defined 
with full weighting (ref. 61), and bilinear interpolation is used for the prolongation operator. 
For full weighting, the restriction operator /£* j is defined by 

^ + 1 ( J +1 — ^ A* j Py ( 1 


(fy+l^T+1 )/ 

= 1 ( n *+l )/ 


where // is a standard averaging operator, and t hus 


and 


f l x R iJ - j( R i+ 1J + 2 R i,j + R i- l,j) 

These operators are adjoint on a uniformly spaced mesh. The operators used in this paper do 
not have the adjoint property. 

Suppose the bilinear-interpolation operator is replaced with a piecewise, constant-prolongation 
operator. This new prolongation operator transfers the same correction to all fine-grid cells that 
comprise a coarse-grid cell. Using the inner product definition of functions, this prolongation 
operator can be shown to be the adjoint of the restrict ion operator of equation (7.3.5). (See ap- 
pendix C.) However, this type of prolongation is not considered an appropriate choice. That is, 
if the Navier-Stokes equations are solved, prolongation does not satisfy the requirement for the 
intergrid transfer operators, which states that the sum of the order nip of the prolongation op- 
erator, and the order m r of the restriction operator must exceed the order 2m of the differential 
operator being considered (ref. 61). With the piecewise, constant prolongation, rnp-j- m r = 2. In 
the case of the bilinear interpolation, nip + m r = 3 > 2 m = 2, and the requirement is satisfied. 
Note that, frequently a restriction operator is chosen that is not the adjoint of the prolongation 
operator. 

7 .3.3. Elements of present method . Section 7.3.2 states that a forcing function is required 
to properly define a coarse-grid problem for the multigrid method. After initialization of the 
coarse-grid solution, the forcing term is constructed as 

Pk^lt+iKk+i-RkOt+iwk+i) 

where Rfc + i is the sum of the residual (_] and forcing function Pfc+\, and 0 < k < A\ If 
k = N — 1, then R\ = /iy. In the case of the multistage time-stepping scheme, the (<y + l)st 
stage becomes 

(9+1) (0) a , "TT 

W L = w l - <y q * i A/j./fy 

where 

R k . = R k .(w { f)+ lf ] 

R k(w ( f) = T + 4M 0) - ad^) 

and the superscripts C and D mean that the discrete operators are associated with the convection 
and physical, viscous terms, respectively. The quantity AD represents the appropriate artificial 
dissipation terms for a given stage. The residuals on G are smoothed with an R-K scheme. 
Information is transferred from one grid to another with a fixed cycle strategy. 

Both V-type and W-type cycles have been considered. Figures 9(a) and 9(b) show the 
structure of these W-type cycles with four and five levels, respectively. The comparable V-type 
cycles consist of the first and last legs of these W-type cycles. Although the W-type cycle 
becomes complex as the number of grids increases, it has a recursive definition. Thus, the 
W-type cycle is essentially as easy to program as the V-type cycle. The work of these cycles is 
as follows for V-type cycle 

4 

Work M(; < - Work FINE; 

and for W-type cycle 

Work \i ( ^ < 2 Work pj\p 

The subscript MG indicates multigrid cycle, and the subscript FINE refers to one time step 
on the finest mesh. The work associated with grid transfer operations has been neglected. At 
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(a) Four- level cycle. 



( b) F i ve-le ve 1 cv cle. 

Figure 9. Structure of multigrid W-type cycle; letter designations are defined as S: solve equations, R: restrictsolution 
and residuals, and P: prolongate corrections. 

a given grid level, additional R-K steps can be performed in both cycles. In particular, the 
application of two R-K steps on and three R-K steps on all successively coarser grids is 

an effective strategy. Multiple coarse-grid time steps reduce the number of cycles necessary to 
reach a prescribed level of convergence (i.e., engineering accuracy, meaning three to four orders 
of magnitude of reduction in the residual). However, the computational time required to realize 
this level is about the same with or without the additional steps. The principal advantage of 
these multiple iterations is the improved smoothing of residuals, which is important for difficult, 
non linear- flow problems. 

Without additional coarse grid sweeps, the W-type cycle generally requires about the same 
amount of computer time for convergence (engineering accuracy) as the V-type cycle. The 
advantage of the W-type cycle is that it provides improved robustness. Therefore, a W-type 
cycle is used in the applications of this paper. 

When solving the Navier-Stokes equations, the viscous terms are computed on each mesh in 
the multigrid process rather than only on the finest mesh (i.e., as in the convective coarse-grid 
correction scheme of Johnson (ref. 62). Computing on each mesh provides improved convergence 
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behavior for low Reynolds number (i.e., 0(1000)) flow cases. For turbulent flows, the viscosity 
associated with Reynolds stresses is evaluated only on the finest grid, and then determined on 
each successively coarser grid by a simple averaging of surrounding finer grid values. Averaging 
is done to obtain a consistent estimate of the eddy viscosity on coarse meshes when an algebraic 
turbulence model is being applied. The artificial dissipation model for the finest grid is replaced 
on coarser grids with a simple, constant-coefficient, second-difference dissipation model. On each 
grid, the boundary conditions are updated at every R-K stage. 

In describing the multigrid method, section 7.3.1 states that on coarse grids approximations 
are constructed for the residual equations at the boundary points (eqs. (7.3.3)). In constructing 
coarse grid approximations, the solution at the boundary points on a coarse grid is driven 
by the residuals for the boundary points on the next finer grid. However, for the present cell- 
centered, finite-volume scheme, such a treatment for the boundary points is not computationally 
convenient (i.e., the boundary points do not lie on the boundaries themselves, but are located 
in auxiliary cells outside of the domain). Instead, the fine-grid boundary conditions discussed in 
section 5 are applied on the coarse grids, and fine-grid accuracy is not maintained at coarse-grid 
boundary points. When transferring coarse-grid corrections to a finer grid, only the changes 
in the solution at the boundary points caused by R-I\ time stepping are used. Although this 
method of treating boundary points can possibly affect the asymptotic convergence rate of the 
multigrid method, it does not change the fine-grid boundary values if the method converges. 

The robustness of the multigrid method is enhanced significantly by smoothing the corrections 
for the fine grid solution. That is, 


M' ( ' ,+1) = W' (H) + AM „H 

where 

AW tot = All j- 4- All r 

The quantity AW j is the solution correction from the finest grid, arid AW C is the resultant 
solution correction front the coarse grids. This smoothing of the corrections reduces the high- 
frequency oscillations introduced by the bilinear interpolation of the coarse- mesh corrections and 
allows convergence of the scheme for a broader range of artificial dissipation coefficients. The 
factored scheme described for implicit residual averaging with constant coefficients ( c ^ = <:,/ = 0.1) 
is used for the smoothing. 

The full multigrid (FMG) method is employed to provide an improved initial solution on the 
finest grid in the multigrid procedure. The FMG method initializes the solution on a. coarser 
grid of the basic sequence of grids, and iterates the solution for a prescribed number of cycles 
using the FAS scheme. The solution is then interpolated to the next finer grid. The process is 
repeated until the finest grid is reached. In this paper, three refinement levels are used for a 
standard mesh density (e.g., 320 by 64 cells). The first and second levels include three and four 
grids, respectively, and 50 cycles are performed on each. There are five grids in the final level. 

8. Turbulence Modeling 

The numerical solution of the instantaneous Navier-Stokes equations for turbulent flows 
requires computing power well beyond what is currently available (ref. 63). To make turbulent 
flow problems tractable using existing computers, a time-averaged form of the Navier-Stokes 
equations must be solved. If the appropriate expansions of Favre variables (ref. 64) are 
substituted for the flow variables in equation (3.1), and the resulting equations are time averaged, 
the n lass- aver aged form of the Navier-Stokes equations is obtained. These equations have the 
same form as their laminar flow counterparts, except that the stress tensor is augmented by the 
Reynolds stress tensor, the heat flux vector is augmented by the heat flux terms associated with 
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turbulence, and additional mean-energy dissipation terms appear (in many cases, these terms 
can be neglected). Closure for this system of time-averaged equations is realized by using the 
eddy-viscosity hypothesis, which states that the Reynolds stress and heat flux terms are related 
to mean flow-field gradients. Moreover, the effective viscosity is obtained by simply adding the 
turbulent viscosity to the molecular viscosity. The Reynolds heat flux terms are approximated 
using the constant, Prandtl-number assumption. Thus, the effective non dimensional transport 
coefficients for diffusion and heat conduction are 


n = HI + Hi 



(8.1a) 

k = kj + ki = 

(p9/ +| 

ca 

(8.1b) 


respectively. The subscript / refers to laminar values, and the subscript t refers to turbulent 
values. The laminar and turbulent Prandtl numbers are 0.72 and 0.9, respectively. 

For aerodynamic computations, the primary requirement for an eddy-viscosity model is to 
provide a good representation of turbulence to allow accurate predictions of mean flow-field 
characteristics. The desire to utilize such capability on a routine basis creates the need for 
turbulence models with a high degree of numerical compatibility. That is, these models must 
demonstrate a favorable interaction with numerical schemes, and must not prevent reliable and 
efficient calculations. This section presents the two turbulent viscosity models applied in this 
paper. Specific modifications of the originally published forms of the models used to improve 
physical modeling and/or numerical compatibility are discussed. 

The basic turbulence model considered is the widely used algebraic model of Baldwin and 
Lomax (ref. 65). This two-layer model defines the nondime nsional turbulent viscosity as 


tH = min [(/i/)-, (fit ) J (8.2) 

where the subscripts i and o denote inner and outer values, respectively. The viscosity in each 
layer is proportional to the product of a length scale and a velocity scale. In the inner part of 
the boundary layer, 

(fit)j = Re p L 2 Q (8.3) 

where Re = Ke/(y/y M), Q is the magnitude of the vorticity vector, L is the length scale given 
by L — KDd , I\ = 0.4 (Von Karman's constant), D represents the Van Driest damping factor, 
and d denotes the distance from the wall. The damping factor D is defined as 


D=1 -exp (-j£) 

d+ = d l ^ P "' iT,) - 

Hw 


( 8.4 ) 


where .4“*" = 26. In this definition of the law-of-the-wall coordinate the original shear stress 
at the wall is replaced with the maximum laminar value. Substituting the maximum laminar 
value prevents the eddy viscosity from vanisliing when the shear stress goes to 0 at a separa- 
tion point. The laminar value eliminates a nonphysical behavior of the turbulence and generally 


62 


removes numerical difficulties. In the outer part of the boundary layer, the turbulent viscosity 
becomes 

(fi t ) 0 = Re C clau C cp p F wake F Kbh (d) (8.5) 

where CQ au = 0.0168 (Olauser's constant), the additional constant C cf) — 1.6, and 

^vvake = m ^ n (^max -^max' ^ wake ^max ) ( 8 -6 ) 

y r max J 

with F max being the maximum value of the function 

F(d) = dQ D (8.7) 


across the layer, and c/ max is the value of d at which F max occurs. The quantity is the 
difference between the magnitudes of the maximum and minimum velocity vectors that occurs 
across the layer. The function /^Kleb(d) represents the Klebanoff intermittency factor, and is 
defined by 


CKIel.(^) 


1+5.5 C\ 


Kleb 


— ) 

diimx J 



( 8 . 8 ) 


Baldwin and Lomax (ref. 65) defined the constant C W ake 1° be 0.25. This value is generally 
unsatisfactory in transonic airfoil flows because it produces oscillatory movement of a shock 
wave. A remedy for this problem is to set C vva k e — 1.0. 

The Baldwin-Lomax (B-L) model just described is also used for wake regions. For wake 
flows, the Van Driest damping factor is set to unity. The B-L model can also be used to 
represent transition to turbulence. However, the specification of a transition location according 
to experiment is generally preferred. 

When implementing the B-L model, care must be exercised when determining the maximum 
of the function F(d), especially for complex flows. Multiple peaks can occur in this function in 
the vicinity of separation. The second peak is chosen in this case. Due to the rapid evolution 
of the numerical solution with the multigrid method, the turbulent viscosity is updated every 
multigrid cycle. Less frequent evaluation can cause either a slowdown or a stall in convergence. 

The B-L turbulence model represents a balance of production and dissipation of turbulence. 
When the boundary layer on a solid surface is subjected to an adverse pressure gradient strong 
enough to cause flow separation, the production and dissipation of turbulence balance break 
down. The inner part of the boundary layer responds immediately to the adverse pressure 
gradient, but the outer boundary layer experiences a delayed reaction. This delayed behavior 
creates a disequilibrium of the two regions. If the size of the separated flow region is large 
enough to alter the surface -pressure distribution, then the history effects cannot be neglected 
in the turbulence model. In general, both the convection and diffusion of turbulence should be 
modeled to accurately predict the turbulent stresses. 

.Johnson and King (ref. 66) proposed a model to account for nonequilibrium effects and used 
the two- layer, algebraic model of Cebeci and Smith (ref. 67) as a foundation for their model. In 
principle, any equilibrium model, such as the B-L model (ref. 65), could be chosen. The basic 
idea of the Johnson- King (J-K) model is to find an appropriate nonequilibrium factor to modify 
the variation of the equilibrium outer-eddy viscosity. The nonequilibrium factor is determined 
so that a transport equation for the maximum shear stress in the boundary layer is satisfied. 
In reference 68, the implementation of a modified version of the J-K nonequilibrium model is 
presented. Reference 69 gives a similar modified form. The elements of these forms of the J-K 
model are described in the remainder of this section. 
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In reference 68, a turbulence reference quantity is defined as 


G 


m — 



(8.9) 


where the index m denotes the maximum of G across the shear layer. This quantity is then 
used to replace the maximum turbulent shear stress divided by density (i.e., the correlation of 
fluctuating velocity components given by appearing in the original J-K model (ref. 66)). 

The advantage of using Gm is that it is invariant with respect to coordinate systems. In the 
formulation of reference 66, the turbulent viscosity is constructed as an exponential blending of 
the inner and outer viscosities. That is, 


Ab 


( a h ) 0 



( 8 . 10 ) 


The inner viscosity is given by 


(p t )j = Re pD 2 K d 



( 8 . 11 ) 


with 


D = 1 — exp 


-d 


\Ti>7 max {p w g m , t w ) 


( 8 . 12 ) 


Ail appropriate valueof .4+ is 17 (rather than the equilibrium value of 26) (ref. 70). The original 
J-K model requires determining the edge of the boundary layer, since the foundation model was 
the Cebeci and Smith model. This requirement is removed in references 68 and 69 bv using the 
B-L model. Moreover, the outer turbulent viscosity is expressed as 


( /b ) 0 & “Tic C C l au C C p pf’ vvake ^Kle b ) 


(8.13) 


where a is the nonequilibrium factor previously mentioned and the other quantities are defined 
the same as for the B-L model. 

Assuming that G is proportional to the turbulent kinetic energy, and introducing a time 
derivative, the ordinary differential equation in reference 66 governing the trajectory of the 
maximum shear stress is replaced with 


dG m 0G m 

+ U rn — 1 " v n 


dG 


di 


dx 


G m 

% + a iT — 
Oy Lni 


G, 


1/2 _ 


(G^) 


1/2 


+ G% 2 V m = 0 


(8.14) 


where u m and v m are the Cartesian velocity components at the location of G rn . The quantity 
(Geq) ni is the equilibrium value of G at the location of G' m , and the length scale L m is defined 
as 


dm 

Lm = 0.4dm < 0.225) (8.15) 

o 

L,„ = 0.09/ (% >0.225) (8.16) 

0 

with 6 being the boundary-layer thickness. An estimate of 6 given in reference 71 is 1.9d max . In 
the original J-K model, the diffusion term V 1U is defined as 


V 


m 


6[0.1-(d m /6)\ 


(8.17) 
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where 


(8.18) 


F(cr) = k^-ll 

The constants a \ and «2 are taken to be 0.25 and 0.5, respectively. With the F(cr) given by 
equation (8.18), there is a singular- like (nonphysical) behavior of the diffusion term at a = 1 . 
In reference 68, F(<t) is expressed as 


F(cr) = max ^0 , cr 1 / 2 — 1 j 


(8.19) 


and thus has a smooth behavior at a = 1. This F(a) makes the diffusion term 0 in regions 
of reverse flow (where <x < 1). The use of equation (8.19) produces greater differences between 
the predicted shock position and the shock position indicated by experiment for a transonic, 
shock-induced, separated airfoil flow w r ith strong nonequilibrium effects. 

if 9 -; 2 is substituted for G m in equation (8.14), the resulting linear equation is given by 


dt 


T u rn 


djhn 

dx 


+ v m 


QOtn 

dy 


+ Sm 



where the source term is 


c — 
Om — 


<n 

2 L„ 


(<Uq) 


-1 


9m 


( 8 . 20 ) 


Equation (8.14) is strictly valid along the curve determined by the maximum shear stress. 
However, equation (8.14) is solved along the solid surface of interest to facilitate t he numerical 
solution method. The misalignment, between these surfaces creates errors in the convection 
terms (ref. 68). To reduce these errors, the velocity components u m and v m are replaced with 
their projections onto the wall boundary. The spatial discretization of the modified equation 
is accomplished by applying the finite-volume technique to the layer of mesh cells adjacent to 
the solid surface. A fourth-difference dissipation term is appended to this semidiscrete equa tion. 
The same five-stage R-K scheme described in section 6 in conjunction with local time stepping 
is used to numerically integrate the equations in time. With implicit treatment of the linear 
source term in equation (8.20), a Courant. number of about 3 can be used. The computation of 
and thus the turbulent viscosity, must be adequately converged to allow convergence of the 
fluid dynamic system of equations. The turbulence model is applied once every time step in the 
solution of the Navier-Stokes equations, and RrK time steps are performed for each update of 
the turbulence field. 

Once the distribution of g m is known, a new variation for the nonequilibrium factor a is 
calculated with the following equation: 


</' +1 = a " 


9 m 


( 8 . 21 ) 


where a n is a at time level n (ref. 70). 


An alternative technique used to solve for g m that is equal to ( — u f v f ),J * is a space 
marching procedure (ref. 66). When applying this procedure, eliminate the time derivative 
of equation (8.20), transform to arbitrary curvilinear coordinates (£.?/). and assume gm is 
independent, of the normal coordinate ?/ (i.e., a £ curve coincides with the transport path of 
(Jm)* Then, obtain 


dg, u _ £ 

< 1 in r tn — “7 

2 


Vm + ~— 


( 8 . 22 ) 
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with U m = £x u m + £y v m- A discrete equation for (<7 m )?+l can be written as 


{ ( $eq ) m [A\ + (flfeq)m. } i (0m)i+ 1 — (^1 + </m)/ T 



(8.23) 


where the index i means evaluated at the previous £ location, 

A t = 

2L m [/ m 

with representing distance along the integration path, and U m denoting the magnitude of the 
velocity vector at the actual location where g m occurs. Again, for convenience, the integration 
path is taken to he coincident with the geometry being considered. This particular approach 
generally seems more robust, and thus is used for all J-K computations. The original argument 
in favor of the time-dependent technique concerned simplicity in extending to three dimensions. 
However, reference 72 indicates that the steady equation for g m can be solved easily with first- 
order, upwind differencing and point-Gauss-Seidel relaxation. Reference 73 provides additional 
discussion on implementation and various forms of the J-K model. 

9. Concluding Remarks 

The elements of a class of explicit multistage time-stepping schemes with centered spatial 
differencing and multigrid are defined and discussed in this report. Additional understanding 
is gained from analysis of a number of components of these schemes. Through this approach, 
the benefit of a local mode analysis in evaluating boundary-point difference stencils for the 
numerical dissipation is demonstrated. The stability of the multistage Runge-Kutta schemes 
Is examined. Hyperbolic and parabolic operator splitting is applied to determine sufficient 
conditions of stability for the Euler and Navier-Stokes (in the absence of convection) equations, 
respectively. The difficulty in rigorously deriving a sufficient condition for the full Navier-Stokes 
equations is discussed. A simple time-step estimate that works well in practice Is given. The basic 
properties of the implicit process of residual smoothing for extending stability are given. Two 
forms of variable coefficients for the residual smoothing procedure are considered. The formulas 
introduced in this report are shown to perform much better than the formulas of reference 20 
for typical meshes used to compute inviscid, airfoil-flow solutions. With these formulas, a new 
set of variable coefficients is constructed that eliminates the general requirement of including a 
diffusion limit in the time-step estimate. The implicit residual smoothing is also used as the basis 
for one of several techniques that are included to enhance the robustness of the basic multigrid 
method. 

Both the equilibrium model of Baldwin and Lomax and the nonequilibrium model of Johnson 
and King are considered for turbulence closure. The implementation of these models, including 
two alternatives for the Johnson-King model, is described in detail. Some modifications to the 
original formulations of the models are made to improve numerical compatibility of the models 
(i.e., make it easier to converge numerical algorithm with the model), and in the case of the 
Johnson-King model, to simplify implementation and improve prediction capability. 


NASA Langley Research Center 
Hampton, VA 2368 1-2199 
January 3, 1997 
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Appendix A 

Equations for Boundary Points 

Consider the elements of the solution vector 

W'=[c v v s] T 


as dependent variables, where c is the speed of sound and $ is the entropy. The Euler equations 
relative to the rotated Cartesian coordinate system (xf, x n ) can be written as 


where 


an d 


at oxt ax n 

A" = A 1 cos 0 -f B' sin 6 
B" — —A' sin 0 + B 1 cos 6 


(Al) 


A' = 


2c( 7 - l)" 1 
0 
0 


c ( 7 - 1)/2 0 0 

u 0 — c ^(7 — 1) _1 

0 u 0 

0 0 u 


B = 


v 

0 

•ir(7 - 1) 
0 


-1 


0 0(7 — 1 )/2 

v 0 


0 

0 


0 v — c 2 (7 — I ) -1 

0 0 v 


In equation (Al), 0 is the angle that the rotated coordinate system makes with the unrotated 
system. Suppose the R.iemann invariants of 1 -D gas dynamics are changed to dependent 
variables. This is done by first assuming t hat the flow is locally homentropic, and by redefining 
the matrices A' and B as the reduced matrices 



u c(y — 1 )/2 0 " 


II 

2c(7-1) _1 u 0 



0 0 v. 



1 

O 

1 

to 


B' = 

o 

o 



_2c( 7 — l ) _1 0 r 

j 


Then, with the matrix 


Q~ 


l 


0 cos 0 

l/v /2 -(7 - 1 ) sin t?/( 2 x/ 2 ) 

.l/v/2 (7 - 1 ) si nt?/( 2 V^) 


— sinfl 

(7 - l)cost?/( 2 v/ 2 ) 
-(7 - l)cos0/(2v/2)_ 


(A 2 ) 
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the reduced form of the solution vector W' can be transformed to a new vector that is a function 
of the Riemann invariants. In addition, the similarity transformation with Q ~ 1 and 


r o 


Q = 


cos 9 


L sin 6 


l/y/2 

— v^ 2 *sin 0 /( 7 — 1 ) 
— \/5*cos0/(7 — 1) 


-1/v^ 

— \Z2sln 0/(y — 1 ) 
-V2 cos 0/(y- 1). 


can be used to diagonalize the reduced form of the matrix B n . Thus, if equation (Al ), with A f 
and B f defined by equations (A2), is premultiplied by Q ~^ , the result is 


<r‘^+ {<r'A"Q) Q~'^+ 


8W 


dW 1 


(A3) 


If Q 1 is considered locally constant, and the variation of W in the tangential direction is taken 
to be negligible, equation (A3) becomes 


dW dW 


= 0 


(A4) 


where A B n is a diagonal matrix of the eigenvalues {</ n . q n + c, q„ — c) of B ” . and W is the vector 
of characteristic variables defined by 


W 


m 


1 7-l p+ 


■R 


1 


iT 


Vi 2 Vi 2 


with 


R + = q n + 

R7 - q„ - 


2 c 

7^1 

2c 

7-1 


qt = u cos 6 + v sin 0 
q n = — u sin# + v cos 6 
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Appendix B 

Development of Residual Smoothing Coefficients 

To obtain insight into an appropriate form for variable smoothing coefficients, consider first 
the approximate factorizat ion scheme 

(/ + ^ m) AW ij = -K;j (B 1 ) 


where AW/ j is the product of the solution vector for the Euler equatioas and the volume Q (as 
determined by a transformation Jacobian), 7 Zjj is the residual vector for the system, and (£, ?;) 
are arbitrary curvilinear coordinates. The operators // and h are standard averaging and central 
difference operators, respectively. Thus, 

= (w j+1/2j - W ; _ |/2J .) 

The transformed, flux-Jacobian matrices are defined jus 


A=l x A + Z„B 

D = 7 ) X A + r/y B ) 

The spectral radii of these matrices are as follows: 

\c 

ff A= a{ ' A) =n 

*2 = 17 


(B2) 


(B3) 


where \ c and A,, are the characteristic speeds defined in equations (4.2.5). If the mat rices A and 
B are approximated as 

,4 — (T~I 

A 

B =<T~ I 

respectively, then equation (B 1 ) can be replaced with 

^ *'/) AW # -j = -Rjj 


A/ 

Q 


1 + k 


(B4) 


when the scalings are taken to be locally constant. Define 

At 

* = in 
A ' = m> 


(B5) 


as the implicit smoothing coefficients for the £ and ?/ directions, respectively. Using equa- 
tions (B5), and taking 
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the smoothing coefficients of equatioas (B5) become 


A| = 


h 1 

+ ^»/ 

h + A '/ J 


(B7) 


Now, consider the case where parabolic implicit smoothing operators are used instead of 
hyperbolic implicit smoothing operators. In particular, consider the implicit finite-volume 
method of Lerat (ref. 52). This scheme includes two stages. The first stage is a physical stage in 
which the change of the solution vector of the Euler equations is evaluated using a Lax -Wen dr off 
scheme. To remove the time step limit of the explicit scheme, a mathematical stage is applied 
in the following integral form: 


// (AW YdQ+Lo—f <r~ fn ■ Y(AW)*j </T = ff A WdQ (B8) 

JJUij 2 Jr. , ,ur. , . - 4 JJ«, , 


'r. i .ur. , 
*+ 


[[ AWffi + w— [ ^ n [n-V(AW](/r= // (AW )*dtt (B9) 

JJq, ; 2 Jr.. i ur. . i B JJi'Lii 

i .)- <- Jr i.i-* l ’J 


t. . iur. . , 

-J+£ ?J_ 7 


where AW is the change in the solution vector, the superscript (*) indicates a provisional value, 
the overbar refers to a value from the explicit physical stage, the quantity ft;j is a mesh cell 
volume, the vector n is a unit normal to the boundary curve I\ and lj is a constant taken to be 
— 1/2. In equations (B8) and (B9), the eigenvalues <r~ and <x~ respectively, are related to the 
spectral radii of equations (B3) as 

a- ft 

= — A . , • 

Y x *i + y >l 

- - *5° 

8 V /x | +!, £ 

Assuming that the quantities inside the integral signs associated with the boundary curves 
are locally constant, and that the curvilinear coordinates £ and r/ are orthogonal, the integral 
equations (B8) and (B9) can be approximated by 


AW* ; 


1 AT 




'dll 2 ' 

_A 

fi 


' + 2 J 


( AW ) * + j j — (AW ) *j 


+ 


1 A t 2 


/<?§r 2 \ 


o 

(AW)?.-(AW)*_, • 


= (AW)-j (BIO) 


i-^rJ 
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< 


[(AW);j + i - (AW)jj] 


AW ij - 


1 A t 2 
4 ftjj 



1 A * 2 
4 Q;j 

where the unknowns are located at the cell centers. If the coefficients 

al r 2 

_A 

Q 


<r~ r 2 

JB 

Q 




= (AW)’ 




(Bll) 




and 

<T ~ | - 

B 

Q 

(which are evaluated at the cell faces) are taken to he locally constant, and the time step is 
defined as equation (136), then the smoothing parameters i3c and ;3 tj depend upon 


Ac 1 

s 

2 

_( Ac + A,/)_ 


A " 1 

2 

.( Af + A,,)_ 

J 


(B12) 


respectively, which are the squares of the smoothing coefficients obtained for the ADI scheme. 

The results from the l-D stability analysis of section 7.2.1, and the understanding of the 
functional dependence of the smoothing coefficients on A^ and A provide a foundation for 
developing /3c and ;3 /; , respectively. To determine formulas for {3 ^ and /3 /; , the 2-D stability of 
a multistage, time-stepping scheme with implicit residual smoothing is examined. Consider the 
2-D, scalar, hyperbolic w r ave equation 


div chv 



(B13) 


Using central difference approximations for the spatial derivatives, a seinidiscrete form for 
equation ( B 1 3 ) is written as 


At 


dw 

IT 



(B14) 


where the Courant. numbers 


At 


At 



By taking the Fourier transform of equation (B14), the following is obtained: 



zw 


n 


where the Fourier symbol z is given by 

£ = —i (Nf sin 6*^ + N tf sin 0 Tf ) (B 16) 

The caret indicates a transformed quantity, and 0 ^ and 0 rf are the Fourier angles for the 
two coordinate directions. If implicit residual smoothing is applied, the Fourier symbol of 
equation (B16) is replaced by 


. N^sinOg + N if sin 0^ 
XtX'i 


where 


\£ = 1 + 20£ (1 - c os 0e ) 

Xf t = 1 + - cos 6 j i ) 

A sufficient condition for stability can be written as 

max \z\ < N* 

for all 0£ and 0 where N* is the Courant number of the unsmoothed scheme. Let 


(bit: 


(B18) 


F — \z I - A r , 


sin 0 c sin#,, 

c + N fj 2- 

\ 0 '/ M\'i 


Then, 


or 


am 


— sin 0e sin 0 ,, 

F < Nc ^ + N fj !L 

X£ Xn 

F < N^f(0c) + N n g(0 J} ) 


e,„ax < A^/max + N,,g max 
Then, a sufficient, condition for stability is given by 

N^/maxT Nij din-dx < N* 

From equation (7.2.7) of section 7.2.1, it follows that 

f - 1 
;,aax _ yrnir 

i 

9 max — — * ~ 

\A + 4 ^/ 

Substituting equation (B19) into equations (B20) yields 

1 1 


N f 


. + N. 


yr+47^ 1 yj ] ~ 4;i, ( 


< A' ’ 


(B19) 


(B20) 


(B21) 
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But 


H r I 

A < = A?7 


A (i 


act. 


AU r 


At, 


= N 


= N 


k 

Xc + A, 


N 
1 +r 







Ac + A, 


-JL.I 

1 + * J 


(B22) 


where » ■,/{ is the ratio of the modified characteristic speeds ( A (/ / Ac ) and is also proportional to 
mesh-cell-aspect ratio. Thus, equation (B21) becomes 


N___\ 1 jV 1 L_<i 

N* 1 + r lf c y / 1 -f 4/^ A* 1 -f \J 1 + 4/t (; 


(B23) 


In the cases of low-aspect-ratio cells (r,^ <C 1) and high-aspect-ratio cells (r ;/ c 1), equa- 

tion (B2.3) can be replaced by 


Ne 1 
— , < 1 

/V* v/1 - 1Ac - 

and 

1±L—1 <1 

N * yr+4^- 

respectively. Thus, write 



(B24) 


Note that these expressions are related to the smootliing coefficients of equations (B12) for the 
implicit method of Lerat. 

The formulas of equations (B24) can also be obtained by substituting the appropriate time 
step estimates into the 1-D smoothing coefficient of equations (7.2.8) in section 7. 2.1. That is, 
the time step of the smoothed scheme is defined as in equation (B6), and the time step of the 
unsmoothed scheme is a 1-D time step. 
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Appendix C 

Multigrid Transfer Operators 

When constructing a multigrid method, appropriate intergrid transfer (i.e., restriction and 
prolongation) operators must be chosen. Often, these operators are selected so that they are 
adjoint operators. Such a choice provides convenience in the analysis of the multigrid scheme 
(i.e., two- level multigrid analysis). In this section, the natural choice for the restriction operator 
of the residual function when a cell -centered, finite-volume scheme is used for discretization is 
considered. Moreover, the restriction process involves simply summing the fine-grid residuals for 
the fine-grid cells that comprise the coarse-grid cell. A piecewise constant prolongation operator 
is shown to be an adjoint operator. 

Consider a 1-D domain S2 = E 3? : 0 < a? < L }. Define a fine grid Gj and a 
coarse grid Gc that, cover the domain such that Gc C Gf. Generate G c by eliminating 
every other mesh point of Gf (delineated by crossed lines in fig. Cl). Let the mesh interval 

(A xj) j = — ^j- 1 / 2 )/ °f Gj be constant. Define h = hf = (Axj)f. Then, the coarse- 

grid mesh interval is h c = 2 h. Let R be the residual function, and let (?;)/, be a correction to 
the fine-grid solution. Assume that the unknowns are stored at the center of a mesh cell. The 
restriction operator for the residual is defined by 

1 2 

'l R '' = rE <'■/«/)< 

C /= 1 

Suppose that the prolongation operator for the coarse-grid correction simply transfers the 
same correction to the fine-grid cells that determine the coarse-grid cell. The operators iff 1 and 
l!^ v are adjoint operators if 

(R>, . 4*2 /,) = [Ui h )*Rh - * 2 /,] = (jfRn . * 2 /, ) (Cl) 


h ^ 


2h 





• 4 • 



1 

1 



1 

1 


Figure Cl. (’ells of two grid levels. 
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where (*, •) denotes an inner product, and the asterisk indicates transpose. To show that these 
operators satisfy equation (Cl), consider the inner product definition for functions. Let the 
index k denote a coarse-grid cell, and let the indices j and j — 1 represent the corresponding 
fine-grid cells. Then 


df/'Rh, v 2h)‘2li = E TTW' 1 + 

k L k k 


(vkhhhk 


and 

(Rh, 4h v 2h)h = E ( R j)h('>khl,hj+ E ( Rj—l )/i ( v k)'2h hj—l 

j even j even 

= [{flj )h hj + {Rj- 1 )// hj- l] ( r k)'2h 

j e veil 

Thus, these operators are shown to be adjoint operators. Note that if the piecewise constant 
prolongation operator is replaced by a linear interpolation operator, the operators are not adjoint. 
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